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The  Xq:  exchange  approximation  is  used  in  self- consistent  APW 
calculations  on  metallic  copper.     The  a  used  in  the  Xa  method  (a=  0.  7225) 
was  chosen  to  yield  zero  pressure  at  the  experimentally  determined  lattice 
spacing.     Self- consistent  APW  calculations  were  carried  out  for  six 
different  lattice  parameters,   using  this  value  of  ot.     The  total  energy 
as  a  function  of  lattice  parameter  resulting  from  these  calculations 
was  used  to  determine  pressure  as  a  function  of  lattice  parameter, 
cohesive  energy,   and  compressibility.     The  cohesive  energy  for  the 
calculation  using  the  experimental  lattice  parameter  and  o;=  0.  7225 
was  found  to  be  0.  286  Ry,  which  is  within  11%  of  the  experimental 
value  of  0.  257  Ry.     The  compressibility  was  calculated  from  two  sets 
of  calculated  pressures  as  a  function  of  lattice  parameter  and  v/as 


i 
\ 

found  to  agree  with  the  experimental  value  to  within  7%  in  one  case  ; 

and  to  within  4%  in  the  other.  j 
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CHAPTER  I 

INTRODUCTION 

1.  1  Purpose 
With  the  development  of  large  fast  computers,   many  methods  of 
determining  the  electronic  structure  of  solids  have  become  practical. 
These  naethods  employ  various  approximations  to  simplify  the  calcula- 
tions, with  a  minimum  loss  of  accuracy.     Within  two  of  these,  namely 
the  nonrelativistic  and  electrostatic  approximations,   the  Hamiltonian 
operator  for  a  crystal  can  be  written  as 

H     =  -T,:kr  v^  -  y  T"  ^^  +  v(x^,  X.)  (1. 1) 

where  the  first  term  is  the  kinetic  energy  term  associated  with  the 
motion  of  the  nuclei  of  masses  M     ,   the  second  is  the  kinetic  energy' 
term  associated  with  the  motion  of  the  electrons  of  masses  m  ,   and 
the  last  term  is  the  potential  energy  term,  which  includes  nuclear- 
nuclear,   electron-nuclear,   and  electron- electron  coulomb  interactions, 
plus  the  electron- electron  exchange  interactions.     In  the  Born- 
Oppenheimer  approximation,   which  is  based  on  the  assumption  that 

since  M     »  m  ,   the  electx^onic  motion  can  be  treated  separately 
M-  o 
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from  the  nuclear  motion,  the  first  term  in  Eq.    (1.  1)  vanishes, 

yieldiag  the  electron  Hamiltonian  as 
»  2 

Hnn--5ZT^      ^^      +       V(X^,X.).  (1.2) 

op       ^-^  2mQ      3  ^      3 

The  Schrodinger  equation  for  the  system  is  then  given  as 

S—-     v^   +  V(X^,x.)    U(X„,x.)  =  E(X    )  U(X^,x.)  (1.3) 

3  2m^      3  ^^      1 J  ^^      D  ^  ^      -^ 

where  the  energy,   E(X  ),   depends  on  the  positions  of  the  nuclei, 

X    ,   only  as  a  parameter.     This  is  the  total  energy'  that  is  considered 

in  this  calculation  as  a  function  of  lattice  parameter. 

Recently,   considerable  effort  has  been  concentrated  on  developing 
an  approximation  to  the  electronic  exchange  interaction  term  in  the 
crystal  Hamiltonian    that  can  be  used  to  determine  the  electronic 
states  of  a  solid  with  reasonable  accuracy.     The  purpose  of  this 
dissertation  is  to  apply  one  such  exchange  approximation  to  the 
evaluation  of  the  energy-bands  of  metallic  copper,  using  the  self- 
consistent  APW  method,   and  to  determine  the  total  energy  as  a  function 
of  lattice  parameter.     These  results  can  then  be  used  to  obtaiji  such 
quantities  as  pressure  as  a  function  of  lattice  parameter,   compressibility 
and  cohesive  energ;^'-. 

The  exchange  approximation  that  is  used  in  this  work  is  known 

JL 

as  the  Xc*:  approximation.     This  method  began  with  the  so-called  P^ 

1  2 

approximation,   originally  proposed  by  Slater.       Later  when  Caspar 

3 

and  Kohn  and  Sham,      using  a  slightly  different  approach,    derived  a 
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1 
p^  exchange  approximation  that  differed  from  that  of  Slater's  original 

result  by  a  constant  niijltiplier  of  2/3,   it  was  indicated  that  a  reasonable 

exchange  approximation  could  be  obtained  by  multiplying  the  Slater 

exchange  by  an  adjustable  parameter,    a.,   thus  the  Xoc  method.     There 

have  been  many  other  exchange  approximations  suggested,   both  local 

and  non-local  in  nature.     However,  the  use  of  several  of  these  on  the 

4 
energy-bands  of  copper  by  Boring  and  Snow    have  failed  to  produce 

results  that  are  significantly  better  than  those  obtained  using  the  less 

complicated  Xa  method. 

This  study  was  conducted  on  copper  for  several  reasons,   the 

first  of  which  is  that  the  author  *     has  done  considerable  work  on  the 

energv'-bands  of  copper'  in  the  past.     There  is  also  a  large  number  of 

other  calculations  on  copper  available  in  the  literature  with  which  to 

7 

compare  and  experimental  results  are  readily  available.     The  author 

has  also  presented  calculations  on  silver,  but  silver,   having  a  greater 
atomic  number  than  copper,   should  be  affected  more  by  relativistic 
effects  which  are  not  considered  in  the  present  work.     Aluminum 
might  have  also  been  considered,   since  aluminum  energy-band 

Q 

calculations  had  also  been  presented  by  the  author,     but  aluminum 

9 
has  been  studied  in  great  detail  by  Ross  and  Johnson.       Another 

reason  for  selecting  copper  over  alum.inum  is  the  large  nimiber  of 

d-like  electrons  in  a  nearly  filled  d-band.     It  has  been  observed  in  past 

work  that  these  d-bands  are  very  sensitive  to  small  changes  in  the 
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exchange  potential  used  in  the  calculation,   whereas  the  broad  sp-band 
as  found  in  aluminum  remains  virtually  unaffected.     This,   of  covirse, 
makes  copper  a  good  test  case  when  considering  various  exchange 
approximations  and  their  effects  on  the  resulting  energy-bands. 

The  method  of  calculation  is  described  in  detail  in  Chapter  II 
with  an  outline  of  the  procedure  given  in  the  next  section.     Chapter  III 
is  devoted  to  the  results  of  the  calcmations  as  applied  to  copper  and 
the  comparison  of  these  results  with  previously  calculated  and  experi- 
mental results.     The  discussions  and  conclusions  are  given  in  Chapter  IV. 

1.2    Outline  of  Calculation 
These  calculations  employ  the  augniented  plane  wave  (APW) 
method       for  determining  the  energy-bands.     In  this  method  the 
"muffin  tin"  approximation  to  the  potential  is  used,  which  has  been 
found  to  give  very  good  resiilts  in  energy-band  calculations  of  this 
type.     In  the  "muffin  tin"  approximation  the  electronic  potential 
energy  function  is  considered  to  be  spherically  symmetric  within 
non- overlapping  spheres  centered  at  each  lattice  site  and  constant  in 
the  region  between  spheres.     DeCicco       has  found  that  the  "non- 
muffin  tin"  terms  should  have  very  little  effect  in  an  energy- band 
calculation  of  an  FCC  metal  like  copper,   and  these  terms  are  not 
considered  in  these  calculations. 
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The  self- consistent  calculation  is  started  with  a  potential  that 
is  generated  by  the  superposition  of  atomic  potentials  centered  at 
each  lattice  site.     This  potential  is  then  used  in  the  APW  calculation 
to  determine  the  energy-bands  and  corresponding  charge  density  of  the 
filled  Bloch  states.     The  energies  and  charge  densities  of  the  core 

states  are  determined  in  a  manner  similar  to  that  used  in  the  Herman  and 

12 
Skillman      atomic  calculations  with  the  appropriate  boimdary  conditions. 

The  total  charge  density,   core  plus  energy-band,   is  then  used  to 
generate  a  new  potential  to  be  used  in  the  next  iteration  of  the  self- 
consistent  process.     This  iterative  process  is  continued  luitil  there  is 
no  significant  difference  between  the  one- electron  energies  obtained  in 
two  successive  iterations. 

This  self- consistent  m^ethod  is  used  for  calculations  at  various 
lattice  parameters  so  that  the  total  energy  as  a  function  of  lattice 
parameter  can  be  determined.     Pressure,   compressibility  and  cohesive 
energy  can  then  be  determined  and  compared  to  experimental  results 
and  the  results  of  previous  calculations. 

These  calculations  are  all  done  for  an  appropriate  value  of  the 
parameter  a  of  the  Xa  method,   the  choice  of  which  will  be  discussed 
in  detail  along  with  the  rest  of  the  method  in  Chapter  II. 


CHAPTER  n 


METHODS  OF  CALCULATION 


2.  1    The  Self- Consistent  APW  Method 
Since  the  introduction  of  the  augmented  plane  wave  (AP"W)  method 

by  Slater       in  1937,  it  has  been  used  for  a  large  number  of  calciolations 

13 
on  solids  with  very  good  results  in  most  cases.     Loucks       and  more 

14 
recently,  Mattheiss,  Wood,  and  Switendick      have  presented  articles 

that  not  only  describe     in  detail  Slater's  original  presentation,  but  also 

cover    the  use  of  group  theory  and  self- consistency.     To  present  the 

APW  method  in  detail  here  would  be  to  rewrite  their  articles.     Since 

this  would  be  a  waste  of  time  and  space,  what  follows  is  a  short  review 

14 
of  the  method  that  parallels  the  article  by  Mattheiss  et  al. 

The  method  is  designed  to  solve  the  on- electron  Schrodinger 

equation  in  a  periodic  potential  V(r)  such  as 

V(r +  T  )    -    V(r)  (2.1) 

where  X     is  any  lattice  vector  of  the  crystal  structure.     The  one- 
electron  equation  is  then  given  by 

H  'lf(iuk)  =U^  +  V(r)|    Mr;k)    =    E(k)    Mr;k)  (2.2) 

vnih  the  energy  B{]^  in  Rydbergs  and  the  distances  expressed  in  Bohr 
radii.     In  the  original  APW  method  as  well  as  in  most  variations  of  it 
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used  recently,  the  potential  functioii  V(r}  is  usually  considered  to  be 
of  the  "nnuiTin  tin"  form.    Although  this  is  not  a  very  strong  restriction 
on  calculations  in  metals  such  as  copper,   DeCicco       has  indicated 
how  one  can  evaluate  the  corrections  associated  with  the  "non-muffin 
tin"  teriTiS  of  the  potential.     The  methods  used  to  construct  the  potentials 
used  in  the  present  work  are  discussed  in  detail  in  Sec.   2.  3. 

The  wave  functions  that  are  solutions  to  Eq.    (2.  2)  are  Bloch 

15 
type  functions       and  can  be  written 

t(£>k)  =  exp  (ik.r)  U(r;;k)  (2.3) 

where  U(r;k)  is  a  periodic  function  of  the  lattice  such  that 

U{r+  T   ;k).  =    U(r;k)  (2.4) 

with  T     equal  to  any  lattice  vector.     It  should  be  noted  that  the  energy 
as  well  as  the  wave  fimction  are  dependent  on  the  selection  of  k ,  the 
wave  vector.     This  vector  corresponds  to  a  vector  in  the  reciprocal 
space  associated  with  the  crystal  structure.     To  solve  the  energy- 
band  problem,   values  of  the  wave  vector  k  must  be  selected,   for  which 
the  solutions  to  Eq.   (2.  2)  are  to  be  determined.     Due  to  the  translational 

symmetry  of  the  recriprocal  space,   one  needs  only  to  consider  vectors 

15 

in  the  first  Bi-illouin  zone.         Also  due  to  the  symmetry  of  the  first 

Brillouin  zone  itself,   one  needs  to  determine  solutions  for  only  those 
wave  vectors  within  1/48  of  the  zone  (in  our  face-centered  cubic  problem), 
as  shown  in  Fig.   1.     Of  course,  this  still  represents  an  infinite 
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number  of  unique  vectors  to  solve  for.     Wlien  one  attempts  to  solve 
an  actual  problem,  the  first  Brillouin  zone  is  divided  into  a  grid  of 
equally  spaced  points,   the  spacing  depending  on  the  accuracy  required 
in  determining  the  resulting  energy- bands.     Again,   one  need  only 
consider  those  points  and  associated  vectors  that  are  contained  in 
1/48  of  the  zone.     Some  of  these  points  and  their  labels  are  also 
indicated  in  Fig.    1.     The  one- electron  Schrodinger  equation,   Eq.   (2.2), 
is  then  solved  for  the  eigenvalues  and  wave  functions  at  each  of 
these  points.     The  resulting  charge  density  associated  with  each  vector 
k  is  then  weighted  by  the  fractional  part  of  an  electron  that  can  be 
accommodated  by  this  state.     This  weight  is  determined  to  be 

W(k)  =  2  •  f,    .48  (2.  5) 

—  K 

where  f,    is  the  fraction  of  the  total  volume  of  the  first  Brillouin 
k 

zone  associated  with  that  point.     This  is  then  miiltiplied  by  48, 
since  only  points  in  1/48  of  the  zone  are  considered  in  the  actual 
calculations  and  the  factor  2  is  introduced  to  accommodate  the  two 
possible  spin  states.     It  should  be  noted  that  if  the  solution  for  the  lowest 
energy  is  foimd  for  each  of  the  points  in  the  first  Brillouin  zone,   this 
represents  a  collection  of  electronic  energy  values,   called  the  first 
band,  that  v/iil  accommodate  tv/o  electrons  of  opposite  spins.     If 
the  potential  has  a  spin  dependent  term,   then  there  are  different 
solutions  for  each  spin  state,        and  the  factor  2  is  left  out  of  the 
equation  for  the  weight,    Eq.    (2.  5). 
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In  the  APW  method,   solutions  to  the  one- electron  Schrodinger 
equation  are  taken  to  be  linear  combinations  of  functions  called  augmented 
plane  waves  (APW's) 

t  fcis)  -  2]v(k.)  $  (r;k.,  E)  (2.  6) 

i 

where  k.  =  k  +  K    with  K    defined  as  a  reciprocal  lattice  vector,   and 
the  v's  are  variational  coefficients.     By  application  of  the  variational 
method,   an  APW  secular  equation  can  be  obtained,   the  roots  of  which 
are  the  approximate  energy  eigem/alues  of  the  state  being  considered; 
solution  of  this  secular  equation  also  yields  the  eigenvectors  v(k^). 

Each  APW  is  defined  as  a  plane  wave  in  the  region  between  "muffin 
tin"  spheres,   called  APW  spheres, 

$  (r^.,  E)  =  exp  (ik|- r)  (2.7) 

and  is  expanded  in  terms  of  spherical  harmonies  inside  the  APW 
spheres.     Within  the  sphere  centered  at  r^, 

CO  Si 

l(rik.,E)=    Y.    E   ^^Jki)U^(Ir-^nl^E)pJ^Us0)exp(imo^(2.8) 

where  the  coefficients  A        (k.)  are  chosen  so  that  the  fimction  is 

j^m  ~i 

continuous  in  value  at  the  APW  sphere  radius  R^.     U^(  ]  r-r^  ]  ;E)  in 
Eq.   (2.  8)  are  the  solutions  to  the  radial  Schrodinger  equation 

._l_-^/,2     ^^>^  ^  /iil±l)  ^.  V  \  u     =EU„  (2.9) 

^2    dr    ^  dr     /       ^  ^2  x)!      H  I 

where  V    is  the  spherically  sym.metric  part  of  the  potential  centered 
n 


about  r   . 
—n 
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Solvirig  for  the  coefficients  A  ^     (k.)  by  matching  the  value  of 

Eq.    (2.  7)  and  Eq.    (2.  8)  at  the  sphere  radius  R   ,   yields  the  APW 

within  the  sphere  centered  at  r 
^  ~n 

f  (rjk.,  E)  =  exp(ik  .^r^)    J2      S  <2^+l>  ^\    (R   -e)       ^i  t"^n  '  '^^ 

£=0    m=  -  i,  X      n 

X   \ll   1^1^^;     pJ"'^cose)P^'"''(cosei)exp[im(f-$.)]  (2.10) 

where  6.  and  §.  are  the  spherical  coordinates  of  k.  with  respect 
to  a  coordinate  system  centered  about  r  . 

The  variational  method,  when  applied  to  Eq.    (2.2),  yields  a  set 
of  N  simultaneous  equations 

N 


y'(H-E)..v(k.)  =  0;       i  =  1,2 

^  ij       1 


N  (2.11) 


Where  N  is  the  number  of  APW's  in  the  expansion  of  the  wave  function 
and  (H-E)..  is  the  matrix  element  between  APW's 

(H-E)..H  <   $(r;k.,E)    lH-EU(r;k.,E)>  .  (2.12) 

Note  that  E  occurs  both  explicitly  and  implicitly.     In  order  to  solve 
for  the  eigenvalues,   one  must  locate  the  energy  E  for  which  the 
determinant  of  the  matrix  of  (H-E)..  vanishes.     This  is  usually  done 
by  evaluating  this  determinant  for  a  set  of  evenly  spaced  values  of  E 
in  the  range  of  interest,   and  locating  the  E's  that  correspond  to  the 
zeros  of  the  determinant  by  linear  interpolation.     The  determinant 
for  this  energ}'^  is  then  evaluated,   and  used  to  locate  the  eigem/alue  more 
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accurately.     This  iterative  procedure  is  continued  until  the  energy  is 
located  as  accurately  as  desired  for  the  calculation  being  considered. 

Group  theory  is  used  to  simplify  the  problem  somewhat.     If 
the  irreducible  representations  of  the  point  group  of  the  wave  vector  k 
are  determined,   they  can  be  used  to  greatly  reduce  the  dimensions 
of  the  determinant  to  be  evaluated,   thus  speeding  up  the  entire  process. 
These  irreducible  representations  also  yield  information  about  the 
transformation  properties  of  the  states  associated  with  the  wave  vector 

k.    The  application  of  group  theory  to  the  APW  method  is  discussed  in 

14 
detail  in  the  v/ork  of  Mattheiss,   Wood,   and  Switendick,        and  will  not 

be  discussed  further  here. 

Once  an  eigenvalue  is  deternriined  to  the  desired  accuracy,   it 
is  used  in  Eq.    (2.11)  to  evaluate  the  eigenvector  v(k.).      With  this 
inforination,   the  spherically  averaged  charge  density  for  the  state 
corresponding  to  k  and  E  can  be  determined  and  weighted  according  to 
Eq.    (2.  5). 

The  eigenvalues  and  corresponding  v/ave  functions  for  the  core 
states,   defined  as  the  lower  energy  atomic-like  electronic  states  that 

can  not  be  determined  by  the  APW  method,   are  determined  in  a 

12 
manner  similar  to  that  used  by  Herman  and  Skillman       in  their 

atomic  calculations.     The  eigenvalues  are  located  by  requiring  that 

the  logarithmic  derivatives  of  the  radial  wave  function  from  outward 

and  inward  integration  of  the  radial  Schrodinger  equation,   Eq.   (2.  9), 

match  at  the  classical  turning  point.     The  boundary  conditions  placed 
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on  the  wave  functions  were 


and 


U^'(R^;E)    =    0;      £  =  0,    2,   4,    .  .  .  (2.13) 

U^     (R^;E)    -    0;      j^  =     1.    3,    5,    .  .  .  (2.14) 


where  the  prime  indicates  the  derivative  of  the  function  and  both  the 
function  and  its  derivative  are  evaluated  at  the  APW  sphere  boundary. 
These  boundary  conditions  were  used  to  preserve  the  translational 
symmetry  and  at  the  same  time  resulted  in  continuity  of  the  value  and 
slope  of  the  wave  function  in  the  direction  of  nearest  neighbors. 

The  total  electronic  charge  density  is  then  determined  as  the 
sum  of  the  weighted  charge  densities  associated  with  the  filled  band 
states  and  the  charge  densities  of  the  core  states.     This  charge  density, 
along  with  that  associated  with  the  nuclei,   is  then  used  to  evaluate  a 
new  potential  and  the  next  step  of  the  self-consistent  procedure  is 
started.     This  iterative  process  is  continued  until  the  difference  in 
energy  associated  with  any  state,  band  or  core,   obtained  in  two 
successive  iterations  is  as  small  as  desired. 

2.  2    The  Xq;  Exchange  Approximation 
In  the  one- electron  Schrodinger  equation  used  in  this  calculation, 
Eq.    (2.  2),   the  spherically  symmetric  potential  can  be  written 

V(r)    =    V^(r)   F  V^(r)  (2.  15) 

where  V   (r)  is  the  one- electron  potential  due  to  the  coulomb  interation 
of  the  electron  under  consideration  with  all  electronic  charges  (including 
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itself)  and  V  ^(r)  is  the  exchange  term.     If  the  Hartree-Fock  form  of 
the  exchange  potential  is  used,   a  different  exchange  potential  is  required 
for  each  electron;  construction  of  these  potentials  involves  sums  and 
integrals  over  all  the  other  wave  fimctions.     In  a  solid  calculation,   with 
the  large  number  of  fiuictions  involved,   a  calculation  using  this  type  of 
exchange  potential  becomes  extremely  difficult. 

What  is  needed  is  an  approximation  to  the  exchange  potential 
which  is  the  same  for  all  electrons  and  easier  to  construct.     Just 
such  a  potential,  based  on  the  average  exchange  potential  in  the 
free  electron  gas,  was  given  by  Slater    as 


[t    '    ^e^>]* 


V^g(r)    =  -  6|^f-    •    P(r)r  (2.16) 

where  p  (r)  is  the  total  electronic  charge  density  at  the  point  r.     If 
p  (r)  is  spherically  symmetric  about  the  origin,   V      (r)  is  also  spherically 
symmetric  and  if  we  define  a  radial  charge  density  p  {r)  by 

p(r)  =   f^  (2.17) 

then  substituting  Eq.    (2.  17)  into  (2.  16),   the  Slater  exchange  potential 
is  given  as 

(2.18) 

2  3 

Caspar     and  later  Kohn  and  Sham    used  a  variational  approach  to 

determine  an  approximation  to  the  exchange  potential  of  the  same  form, 
but  found  it  to  be  2/3  as  large  as  the  Slater  exchange  approximation. 
This  indicated  that  possibly  a  reasonable  approximation  to  the  exchange 
potential  would  be  wJiat  is  now  called  the  Xo:  exchange  approximation, 
given  by 
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where  a  (a  =  2/3  in  the  Kohn-Sham  case)  is  an  adjustable  parameter 
to  be  selected  to  fit  the  problem  under  consideration.     Since  then, 
several  calculations  have  been  presented  on  both  atoms  and  solids  using 
this  parametrized  one- electron  exchange  potential,  with  reasonable 
success. 

There  is  still  the  question  of  the  choice  of  the  parameter  a.  for 

17 
a  particular  calculation.     In  a  recent  paper,  Schwarz       describes  three 

meftiods  for  determining  a.  in  atomic  calculations.     The  first  was 

1  Q 

suggested  by  Lindgren,        in  which  the  a  is  chosen  so  that  the  Hartree- 

Fock  total  energy,   determined  from  the  Xa:  orbitals,   is  a  minimum. 

This  a  ,   called  a.     .    ,   has  been  evaluated  for  several  elements 
mm 

3  9  20 

by  Kmetko       and  Wood  .        The  second  method,   suggested  by  Berrondo 

21 

and  Goscinski,        is  based  on  the  quantum- mechanical  virial  theorem 

and  recommends  that  the  a   ^  should  be  chosen  so  that  the  virial 

vt 

coefficient  T],   given  by 

Tl  =    -  <P.  E.  > /2<K.E.  >  (2.20) 

where  <P.  E.  >  and  <K.  E.  >  are  the  average  potential  and  kinetic 
energies,   respectively,   calculated  from  the  Hartree-Fock  expression 
using  'Xci  orbitals,   be  equal  to  1.     The  third  method  chooses  Q,   called 
a       ,   so  that  the  statistical  total  energy,   defined  for  a  crystal  in 

Sec  2.  4,   of  the  Xc  method  equals  the  Hartree-Fock  average  energy 

22-24 
of  configuration  for  the  ground  state  of  the  atoin.     Schwarz  has 
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determined  the  values  of  C^   ,  and  cy^^,-,  for  the  first  41  elements  and 

vt  HF 

has  found  that  they  differ  by  less  than  0.  0008  in  every  case  except 

25 
hydrogen.     Slater  and  Wood  '  have  recommended  using  aTxr-.  i^^  ^ 

HF 

solid  calculation.     Energy  band  calculations  similar  to  the  one  presented 

9  R    9  7 

here,   using  values  of  a  given  by  Schwarz  have  been  done  by  Averill     ' 
and  Hattox. 

In  the  present  calculation,    a  is  chosen,   that  is  similar  to  a    .  > 

VT 

such  that  the  virial  coefficient  equals  unity  in  the  crystal  at  the  experi- 
mentally determined  lattice  parameter  appropriate  for  0°K.     This  is 
equivalent  to  requiring  that  the  minimum  of  the  total  energy  vs.  lattice 
spacing  curve  fall  at  the  experimental  value  of  the  lattice  parameter 
(A  )  as  shovm  in  Fig.   2;  2.      Since  the  APW  method  already  requires  the 
lattice  parameter  as  an  input  parameter,   this  method  of  determining  a 
adds  no  new  calculational  parameters  (either  from  experiment  or  from 
atomic  calculations).     There  are  two  ways  to  find  this  ci  •     One  is  to 
determine  the  movement  of  the  minimum  of  the  total  energy  vs.  lattice 
spacing  curve  as  a  function  of  a  and  then  to  use  this  information  to 
determine  the  desired  value.     This,   of  course,   requires  a  reasonably 
large  number  of  self- consistent  APW  calculations  for  the  various  values 
of  c.  and  the  lattice  parameter. 

The  other  way  involves  the  virial  theorem,   given  by 


<K.E.  >  -  -  1/2 


.       1    1 


(2.21) 


where  <K.   E.  >  is  the  average  kinetic  energy,  x.  are  the  coordinates 

of  tlie  system  and  F^.  ^g  corresponding  component  of  force,   arising  botii 

29 
from  forces  within  the  system  and  external  forces.     Averill       has  shov/n 
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Cohesive  Energy 

1 


Lattice       Parameter 


Figure  2.  2.     Qualitative  graph  of  the  ground  state  total  energy 
of  a  crystal  as  a  function  of  the  lattice  parameter. 
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that  for  the  Xa  method  in  the  "muffin  tin"  approximaticn ,   the  forces 
within  the  system  result  in  a  term  given  as 

<  P.  E.  >  (2.22) 


fc^-v] 


av 
1 

where  <P.  E.  >  is  the  average  potential  energy,   and  the  external  forces 

result  in  a  term  representing  the  volume  times  the  pressure  external 

to  the  system  which  is  required  to  "hold"  it  in  the  assumed  configuration. 

The  pressure,   P,   is  given  by 

P=-^       <P.  E.  >  +  2  <K.  E.  >|  (2,23) 

where  V  is  the  volume  and  P  =  0  corresponds  to  the  virial  coefficient 

equal  to  1.     This  pressure  is  then  determined  as  a  function  of  a.  and 

that  0!  selected  that  yields  zero  pressure  at  the  experimental  lattice 

parameter.     This  does  not  require  quite  so  many  calculations,  but  the 

calculations  are  much  more  difficult  to  carry  to  convergence;  the 

requirement  is  now  that  the  calculated  pressure  be  stable  upon  proceeding 

from  one  iteration  to  the  next.     Obviously,   one  method  is  about  as 

involved  as  the  other.     However,   the  latter  method  was  used  in  the 

calculation  presented  here. 

2.  3    Potentials 
The  crystal  potential  used  in  these  self- consistent  APW 
calculations  is  of  the  "muffin  tin"  type,   as  mentioned  in  Sec.    1.  2; 
that  is,   it  is  spherically  symmetric  within  the  non- overlapping  spheres 
surrounding  each  atom  and  constant  between  spheres.     The  radii  of 
these  spheres,   called  APW  spheres,   are  equal  to  one  half  the  first 
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nearest  neighbor  distance.     This  insures  that  the  spheres  will  contain 
the  largest  fraction  of  the  Wigner-Seitz  cell  volume  associated 
with  each  atom,  without  overlapping.     In  practice,   it  is  found  that  the 
larger  this  fraction,  the  better  the  "muffin  tin"  potential  approximates 
the  exact  one. 

The  coulomb  part  of  the  starting  potential  used  in  the  self- 
consistent  procedure  is  generated  by  a  superposition  of  the  coulomb 
part  of  the  atomic  potentials  of  the  first  few  nearest  neighbors, 
spherically  averaged  about  the  central  atom.     This  is  achieved  by 
expanding  the  spherically  symmetric  atomic  coulomb  potentials  of 

the  neighboring  atoms  about  the  center  atom  by  use  of  the  alpha- 

30  31 

function  expansion  suggested  by  Lowdin,        as  described  by  Mattheiss. 

In  metals,   including  contributions  from  the  first  four  nearest  neighbors, 

usually  gives  the  required  convergence  of  this  expansion.     The  coulomb 

part  of  the  constant  potential  between  APW  spheres  is  set  equal  to  the 

average  of  the  superposed  coulomb  potential  in  the  region  between 

the  sphere  and  the  Wigner-Seitz  cell  boundary. 

The  exchange  portion  of  this  starting  potential  is  determined  by 

the  Xa  approximation,   using  a  crystal  charge  density  in  which  the 

contribution  of  neighboring  atoms  have  been  included  by  use  of  the 

same  expansion  used  to  expand  the  superposed  covilomb  potential. 

The  excha.nge  term  in  the  constant  potential  is  determined  by 

V    -eaf-f-   p  ]*  (2.24) 

ox  L8rr  cj 
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where  p     is  the  constant  part  of  a  ''muffin  tin"  charge  density  and  is 

determined  by  requiring  each  Wigner-Seitz  cell  to  be  electrically 

neutral. 

To  determine  the  potential  for  the  second  and  following  iterations, 

the  charge  density  resulting  from  the  previous  iteration  is  forced  into 

a  "muffin  tin''  form,,  as  mentioned  above.     The  coxilomb  part  of  the 

potential  inside  a  sphere  is  evaluated  using  the  spherically  symmetric 

part  of  the  charge  density  in  an  integral  solution  to  the  Poisson  equation 

given  as  ^^ 

^c^''^  "  "  "^  "^r    J    P(r')dr'-2J       -^^  dr'  +  C  (2.25) 

0  r 

where  Rg  is  the  APW  sphere  radius  and  C  is  a  constant  of  integration, 

needed  only  if  the  potential  function  is  to  be  located  with  respect  to 

the  absolute  zero  of  energy.     The  total  potential  function  contains  a 

finite  discontinuity  at  Rg,  the  coulomb  part  of  which  is  determined 

32 

by  use  of  the  Ev/ald  potential  method  as  described  by  Slater  and  DeCicco.  " 

Solution  of  this  problem  insures  that  account  is  taken  of  the  effects  of 
all  nuclei  and  "muffin  tin"  charge  densities  located  exterior  to  the  cell 
for  which  V     is  being  computed. 

The  exchange  term  is  determined,  as  in  the  superposed  case, 
by  the  Xa  approximation.     The  discontinuity  in  this  term  at  Rg  is 
defined  as  the  difference  between  the  exchange  corresponding  to  the 
spherically  symmetric  charge  densit>^  within  the  spheres  evaluated  at 
RS  and  the  exchange  corresponding  to  the  tmiforni  charge  distribution 
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between  spheres. 

The  total  potential  function  is  then  given  as  the  sum  of  the  cotilonib 
and  exchange  terms, 

V'(r)    =    V   (r)    +  V    (r)  (2.26) 

ex 

with  a  total  discontinuity  at  Rg  equal  to  the  sum  of  the  coulomb  and 

exchange  discontinuities.     The  potential  V.  used  in  the  APW 

calculations,   for  convenience,   is  shifted  by  a  constant  V     so  that  the 

constant  potential  between  spheres  is  zero,   and  is  given  by 

^APW^^^    =  V'(r)^V^'  (2.27) 

where  the  primes  are  used  to  indicate  that  these  potentials  are  not 
located  with  respect  to  the  absolute  zero  of  energ}^  at  infinity  as  is 
the  case  of  the  atomic  potential  function,   and  the  case  of  the  super- 
posed potential  used  as  a  starting  potential. 

DeCicco       has  suggested  a  method  to  correct  this  potential 
such  that  it  v/ill  be  in  terms  of  the  absolute  values  of  energy.     This 
method  is  based  on  the  argument  that  the  change  in  the  potential  at 
any  point,   in  going  from  the  intial  to  a  final  interation,   is  completely 
determined  by  the  change  in  the  charge  distributions.     The  change  in 
the  coulomb  part  of  the  potential  at  any  point  is  given  as 

6V(r)  =  2    I   f  ^^^,'i  d^r'  (2.28) 


where  the  integration  is  to  be  performed  over  the  entire  solid.     However, 
with  the  proper  choice  of  cells  over  which  the  integral  is  evaluated,   it 
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can  be  made  to  converge  rapidly.    A  further  simplification  is  achieved 
by  evaluating  the  change  in  the  coulomb  potential  at  a  point  between 
spheres  where  only  the  Ewald  charge  distribution  (assuming  6  p  is  of 
"muffin  tin"  form)  contributes  to  the  constant  potential  and  need  be 
considered.     This  will  yield  the  change  in  the  constant  part  of  the  coulomb 
potential  between  the  initieil  and  final  iteration,  which  is  just  the 
correction  needed  to  locate  the  final  potential  with  respect  to  the 
absolute  zero  of  energy. 

It  can  be  shown  that  if  the  change  in  coulomb  potential  is  evaluated 
at  r  =  (lll)A   /4  in  the  FCC  structure,   the  Ewald  method  of  solving  the 
potential  problem  can  be  used  to  yield 

^'^0  =  ^Oe  -"  5  v((lll)A^/4)  -  V^    ((lll)-\/4)  (2.  29) 

where  V       is  the  change  in  the  average  Ewald  potential  between  spheres, 

V  '  ((lll)A„/4)  is  the  change  of  the  Ewald  potential  at  r  =  (lll)A    /4, 
e  0  u 

and  6V  ((lll)A   /4)  is  the  change  in  the  coiolomb  potential  at  the  point 

{111)A   /4  in  going  from  the  initial  to  the  final  iteration.     The  values 

of  V        and  V     ((lll)A^/4)  for  an  FCC  structure  are  determined  by 
Oe  e  0 

32 
Slater  and  DeCicco       to  be  given  by 

V    '  =   -^-P      (0.8163)  (2.30) 

Oe  Aq 


and 


V^'  =   ((IIDAqM)  =-2^      (0.8019)  (2.31) 


0 

where  A    is  the  lattice  constant  in  Bohr  radii  and  6Q  is  evaluated  by 
subtracting  the  Ewald  charge  for  the  initial  iteration  from  that  of  the 


•24- 


final  itei-ation.     DeCicco       has  also  evaluated  6V  ((lll)A   /4)  for 
an  FCC  structure  and  found  it  to  be  given  by 

6v((lll)A^/4J     =      ^-^    (0.2796).  (2.32) 

Substituting  these  values  into  Eq.    (2.  29)  yields 

^V_     ="T^     (0.8163  +  0.2796-0.8019) 
Oe       A^ 

or 

^^Oe  ""^A^    (0.2940).  (2.33) 

The  constant  coulomb  potential  V       for  the  final  iteration,  located 
with  respect  to  the  absolute  zero  of  energy,   is  then  given  by 


where  V       (SFA)  is  the  coiolomb  part  of  the  constant  potential  for  the 
superposed  free  atom  case.     The  total  constant  potential  is  then  given 


V     =  V       +  V  (2    3  5) 

0  Oc         Ox  K^.-iO) 

where  V    ^  is  the  exchange  part  of  the  constant  potential. 


2.4    Total  Energy  and  Cohesive  Energy 

If  the  Xa  approxim.ation  to  the  exchange  potential  is  used,   the 

33 

expression  for  the  total  energy  of  a  crystal,   as  described  by  Rudge, 

may  be  written  as 
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E  =  V 


E/d^r^  U:^(r^) 


-V^    +      V      (rj+^  V      (rj 
n.e  •^l        2      ee  ~1 


2    2  ex  ^1 


1  ~1 


(2.36) 


where 
V 


~J2  15 


z  z. 

a b 


a,b 
a/b 


R  -R 
a   -^b 


(Nuclear- Nuclear  Interaction) 


2 

-  V  =  (Electron  Kinetic  Energy) 

2Z 

V      (rj  =/  ^  I  -n  ■         r  =  (Nuclear- Electron  Coulomb  Potential) 
ne  ~1      ^— '    R   -r    1 
a    I  ~a  ~1  ' 


V      (r  J  =   /d  1-       , 
ee  ~1       J      ~2    Ir^  -  r^ 


V.      (rJ  =  -67^  P(rj 
ex  "^1  8n         -^1 


(Electron- Electron  Coulomb  Potential) 


(Average  Free  Electron  Exchange) 


and 


p(rj  =y2uT  (r   >  U.(r  )         =  (Total  Electron  Density) 


The  one- electron  Schrodinger  equation  can  be  used  to  eliminate  the 

2 

-V     term  and  the  expression  rewritten  as 


E  = 


/O  0  -, 

d^r    fv      (rJ  +  V      (r  J  +  aV       (rj     p(r) 
~1  L    ne  ~1  ee~l  ex   ~1  J      ~1 

-       L 

fi   fd^^i  [v      (r  J  +  V^   ( r   )  +1  a  V      (r  ■)]    p  (r   ) 
I  2  J      ~1  L    ne  ~1  ee     i        ^  ex  ~i  J       ~i  _ 

|-J  /d^^V      (rj   p(r.)  +  V      1 
\_2J      ~1    ne  ~1        ~1  nn  I 


(2.37) 


where  the  superscript  0  on  Y       and  V      indicates  that  these  potential 
^  ^  ee  ex 

terms  are  evaluated  for  the  charge  density  from  the  previous  iteration. 
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whereas    the  other  potentials  terms  are  derived  for  the  charge  density- 
resulting  from  the  final  iteration.     With  Eq.    (2.  37)  grouped  this  way, 

the  first  term  is  the  average  kinetic  energy  of  the  electrons  (KE)       and 

av 

the  last  two  represent  the  average  potential  energy  (PE)      .     It  has 

av 

been  shown  by  DeCicco       tliat  the  last  term  can  be  written  as 

where  Va  o  is  the  electrostatic  potential  at  the  ath  nucleus  due  to 
all  of  the  other  charge  in  the  crystal. 

Because  of  the  translational  symmetry  of  the  crystal  structure, 
this  expression  can  be  written  in  terms  of  N,   the  number  of  total  cells 
in  the  crystal,   times  the  total  energy-  per  Wigner-Seitz  cell.     The 
expression  for  the  total  energy  per  cell,   in  the  cass  of  one  atom  per 
cell,   can  be  written 


2^e'-    /   V(r)p(r)dr     +     ^        \v     (r)  +  V      (r) 
^1    '      -/cell  J  ^Jcell''^  ^^ 

+  |aV^^(r)jp(r)dr+f     V^_ 


cell 

(2.39) 


where  p  (r)  is  the  radial  electronic  charge  density  and  the  primes  on 
C^'  and  V'(r)  indicate  that  both  are  evaluated  with  respect  to  zero 
constant  potential  between  spheres.     The  potentials  in  the  second  term 
can  be  written 


V      (r)  =    -  -2^ 
ne  r 


I      p(r')dr'+2J     ""  "^^^^  dr'  +  C 


V      (r)  =  - 
ee  r    ,  . 

0  -^  r 
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V      (r)  =  -  ef-^    p(r)l^ 
ex  L8tt  J 


and 


Vq  =  -2   t^s  P^^f'^      dr'    +C  (2.40) 


: 


0 

where  C  is  the  undetermined  constant  of  integration  that  represents 
the  effects  of  all  the  charge  outside  of  the  APW  sphere.     The  V 
term  can  be  written 


fvo^-zf^-^*.    .f  c 


1o         '•■ 


while  the  term  involving  C  in  the  preceding  integral  is  given  by 


p(r)  dr    =      --f-    C.  (2.41) 

'cell    .  ^ 


Therefore,   the  terms  depending  on  the  constant  value  of  the  potentials 
cancel  and  the  total  energy  is  independent  of  the  absolute  zero  of  energy 

for  the  potential  functions. 
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The  cohesive  energy  (E  )  is  defined  as 

E     =  E^.     -     E/N  (2.42) 

c  FA 

where  E     .    is  the  ground  state  total  energy  of  the  free  atom.     This  is 

shown  graphically  in  Fig.   2.  2  of  Sec.   2.  2  above. 

2.  5    Pressure  and  Compressibility 
There  are  two  ways  to  determine  the  pressure  in  a  crystal. 
The  first  method  involves  the  use  of  the  virial  theorm  associated  with 
the  calculation  for  a  metal,   in  which  both  the  Xa  and  "muffin  tin" 
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approximations  are  used,   as  decribed  in  Sec.   2.  2.     With  this  method 
the  pressure  is  given  in  Eq.    (2.  23)  as 


p  =  ~      2  (KE)^     +(PE) 
3V  av  ay. 


.]. 


If  the  (KE)       and  (PE)       are  expressed  inRydbergs  (Ry),   and  V  is 

av  av  ^  J  a         J  ' 

expressed  in  cubic  atomic  units  (cau),  the  pressure  in  kilobars  (Kb) 
is  given  by 


P  =    (4.8959  X  10^)   — 


2  (KE)       +  (PE)  .  (2.43) 

av  av 


The  second  method  of  determining  the  pressure  is  to  generate  the 
curve  for  energy  vs.   volume,  by  evaluating  the  energy  for  several 
values  of  the  lattice  constant.     The  pressure  is  then 

•    P=-fr  (2.44) 

which  can  be  evaluated  for  any  point  on  the  cuive  by  numerical  methods. 
As  in  the  first  method,   the  pressure  in  kilobars  (kb)  is  given  by 

P  =(-1.4688  X  10^)-^  (2.45) 

dV 

with  the  energy  and  volume  expressed  in  atomic  lonits. 

The  compressibility  (K)  given  as 

K  =  --^    ^-  (2.46) 

V      dP  \^.^"/ 

is  equal  to  the  reciprocal  of  the  bulk  modulus,   given  by 

B=Tr  =  -^^  •  <2-«' 

Substituting  for  P  in  this  equation  from  Eq.    (2.44),  yields 

£e^  (2.48) 

B  -  V  ^^,2      . 
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Therefore,   the  compressibility  can  be  evaluated  at  any  value  of  the 
lattice  constant  by  numerically  evaluating  the  second  derivative  of  the 
energy  vs.   volume  cu.rve  generated  for  use  in  the  second  method  of 
determ.ining  pressxire,   described  above. 


CHAPTER  in 

RESULTS  OF  THE  CALCULATION  FOR  COPPER 

3.  1    Determining  oi  of  the  Xoc  Method 
The  a.  used  in  this  calculation  was  determined  to  be  the  one  that 
v/ould  yield  a  virial  coefficient  equal  to  unity,   as  described  in  Sec.   2.  2. 
This  was  done  by  evaluating  the  pressure  for  a  series  of  self- consistent 
APW  calculations,   each  for  a  different  value  of  a.     The  pressure  as  a 
function  of  o.  from  these  calculations  was  then  plotted,   and  the  a. 
determined  that  would  yield  zero  pressure  at  the  experimentally  observed 
lattice  constant  (A^.  The  pressure  for  each  calculation  was  determined  by 

(3.1) 

which  comes  from  the  virial  theorem  as  described  in  Sec.   2.  5. 

Each  self- consistent  calculation  was  done  using  a  500  point 
linear  radial  mesh  (for  integration  of  the  radial  Schrodinger  equation) 
and  for  2048  points  in  the  first  Brillouin  zone.     The  self- consistency 
criterion  used  in  these  calculations  was  that  the  resulting  pressure 
vary  by  less  than  1  Kb  between  two  successive  iterations.     This 
required  from  8  to  10  iterations  in  most  cases,   with  each  iteration 
takjjig  about  5  minutes  of  computer  time  on  tb.e  CDC  7600.     It  is  of 
interest  to  note  that  the  variation  in  total  energy  between  two 
successive  iterations  was  about  10    '  Ry  when  the  pressure  criterion 
v/as  satisfied. 
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P    =    -^ 

3V 


2  (K.  E.  )       +  (P.  E.  ) 

av  'av 
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The  resulting  pressures  for  four  such  calculations  are  given 
in  Table  3.  1,   along  with  the  kinetic,   potential,   and  total  energies. 
Pressure  as  a  function  of  o:  is  plotted  in  Fig.   3.  1,   while  total  energy 
as  a  function  of  c*:,   for  the  same  four  calculations,   is  given  in  Fig.    3.  2. 
The  total  energy  is  very  nearly  a  linear  fimction  of  c,   as  is  the  case 
in  atomic  X©  calculations,  while  the  pressure  is  not.     However,   a 
smooth  curve  can  be  drawn  through  the  average  of  tliese  data  points  on 
the  pressure  vs.    a  plot  in  Fig.    3.  1;  this  predicts  an  a.  of  about  0.  7225 
for  zero  pressure.     This  was  the  value  of  a.  used  to    evaluate  the  total 
energy  as  a  function  of  lattice  parameter  and  the  resulting  pressures 
and  compressibility. 

3.  2    Total  Energy  and  Cohesive  Energy  Results 
Calculations  were  m^ade  for  six  values  of  the  lattice  parameter  (A) 
of  FCC  copper,  each  with  a=  0.  722  5.     These  were  the  lattice  parameter 
corresponding  to  the  experimentally  determined  Wigner-Seitz  cell 
volume  (V  )  and  lattice  parameters  corresponding  to  10%  and  20% 
reductions  of  this  volume  and  10%,   20%,   and  60%  expansions  of  this 
volume.     Each  calciolation,   with  the  exception  of  the  60%  expansion, 
was  converged  until  the  pressure,   as  determined  using  the  virial 
theorem  (Eq.    2.  43),   varied  by  less  than  1  Kb  between  two  successive 
iterations.     This  usually  took  from  10  to  15  iterations,   at  five  minutes 
an  iteration  on  a  CDC  7600  computer.     As  in  the  calculations  discussed 
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TABLE  3.  1. 


The  pressure,  potential  energy,  kinetic  energj^   and  total  energy  for 
four  self-consistent  APW  calculations  on  copper,   with  different  values 

of  (X' 


Pressiare 

Potential  Energy 

Kinetic  Energy 

Total  Energy 

& 

(Kb) 

(Ry) 

(Ry) 

(Ry) 

1 

0.70000 

33.70 

-6553.880469 

3276.967354 

-3276.913115' 

i 

0.70635 

24.35 

-6556.226508 

3278. 132847 

I 
-3278.093661 

0.72000 

3.C4 

.  -6561.271765 

3280.638812 

-3280.632953 

0.77000 

-60.34 

-6579.  806760 

3289.854829 

-3289.951931; 

Figure  3.  1.     Pressure  vs.   a  of  the  Xa  method,   for  copper. 


-34- 


60 


40 


20 


-20 


-40 


-60 


-80 


_L_ __J_ 

.70  .75 

a  of   X„  Method 


Figure     3.1 


Figure  3.  2.     The  total  energy  vs.    a.  of  the  Xo:  method  for  copper 
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in  the  preceding  section,   when  this  convergence  criterion  was  met, 

-  6 

the  change  in  total  energy  was  about  10       Ry  between  two  successive 

iterations.     The  calculation  corresponding  to  60%  volume  expansion 

-4 
was  converged  so  that  the  total  energy  varied  by  less  than  10       Ry 

between  two  successive  iterations.     Each  calculation  was  done  for 

2048  points  in  the  first  Brillouin  zone  and  used  a  500  point  linear 

radial  mesh. 

The  total  energies  resulting  from  these  calculations  are  given 
in  Table  3.  2.     Also  given  are  the  volumes  in  cubic  atomic  units  (cau) 
and  the  lattice  parameters  in  atomic  units  (au).     These  data  are 
plotted  in  Fig.    3.  3,   showing  the  total  energy  as  a  function  of  lattice 
parameter.     This  curve  shows  a  definite  minimum  near  the  experimentally 
determined  lattice  parameter  (A-_),   as  predicted  by  the  results  given  in 
the  preceding  section  (Sec.   3.  1). 

The  cohesive  energy  was  determined  as  the  difference  in  the 
total  energy  at  the  minimum  of  the  energy  vs.   lattice  parameter 
curve  in  Fig.    3.  3  and  the  statistical  total  energy  for  the  atomic 

calculation  using  the  same  value  of  Q,   as  described  in  Sec.   2.4.     The 

2 
atomic  calculation  for  the     S  ground  state  was  a  non- spin-polarized 

12 

Herman- Skillman       calculation  using  the  Xa  exchange  approximation. 

These  total  energies  are  approximately  -3281.  098  Ry  and  -3280.  813  Ry 

for  the  metallic  and  atomic  calculations,   respectively.     This  gives  a 

34 
cohesive  energy  of  0.  286  Ry,   compared  to  the  experimental  value       of 
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TABLE  3.  2 


Total  enei-gy  for  six  values  of  the  lattice  parameter  for  copper  (o;=0.  7225). 


v/v^ 

Volume  (V) 
(cau) 

Lattice  Parameter  (A] 
(au) 

Total  Energy 

(Ry) 

0.80 

63. 1256584 

6. 3205563 

-3281.074352 

0.90 

71.0163639 

6. 5736434 

-3281.093786 

l.'OO 

78.9070739 

6.8086129 

-3281.098246 

1.10 

86. 7977809 

7.0283957 

-3281.094903 

1.20 

94. 6884870 

7. 2352308 

-3281.087333 

1.60 

126. 2513250 

7.9634021 

-3281.044053 

03 

(Free 
Atom.) 

00 

00 

-3280.812553 
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0.  257  Ry.     The  fractional  difference  between  this  calculated  value 
and  the  experimental  value  of  cohesive  energy  is  0.  112  or  about  11%. 
This  is  fairly  good  agreement  if  one  considers     that  the  calculated 
value  of  cohesive  energy  is  about  0.  00001  of  the  total  energies  used 
to  determine  it. 

3.  3    Pressure  and  Compressibility 
The  total  energy  as  a  function  of  volume  is  shov/n  in  Fig.    3.  4 
for  the  first  five  calculations  given  in  Table  3.  2  of  the  preceding  section. 
This  covers  only  a  small  part  of  the  entire  curve  near  the  minimum 
showing  how  the  slope  and  curvature  vary  in  that  region. 

The  pressures  at  each  of  these  five  values  of  the  volume  for 
these  calculations  Avere  determined  using  two  diffex-ent  metriods  and 
are  given  in  Table  3.  3.     The  pressures  given  in  the  third  column  of 
this  table  were  determined  using  Eq.    (2.  43)  derived  from  the  virial 
theorem  as  described  in  Sec.   2.  5.     The  pressures  given  in  the  fourth 
column  were  determined  as  the  nega.tive  of  the  first  derivative  of  the 
total  energy  as  a  fiuiction  of  volume  as  given  in  Eq.    (2.44)  (evaluated 
for  each  value  of  the  volunae^   also  described  in  Sec.   2.  5).     Both 
sets  of  pressures  are  sho\vn  in  Fig.    3.  5,  plotted  as  a  function  of 
volume.     The  agreement  between  the  pressures  determined  by  these 
two  methods  is  obviously  quite  good. 

The  bulk  modulus  was  e\'aluated  using  both  curves  shown  in 
Fig.    3.  5,  and  is  determined  as  the  volume  times  the  negative  of  the 
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TABLE  3.3 

The  pressure  for  five  values  of  the  cell  volume  of  FCC  copper,   evaluated 
using  the  virial  theorem  (vt)  and  the  slope  of  the  energy  vs.   volume  curve 
at  Vq  (    E  vs.  V).    (Vq  is  the  volume  corresponding  to  the  experimental 
lattice  parameter. ) 


V/Vq 

Volume  (V) 
(cau) 

Pressure  (vt) 
(Kb) 

Pressure  (E  vs. 
(Kb) 

V) 

o.'so 

63. 1256584 

559.  9 

562.3 

0.  90 

71.0163639 

200.9 

194.6 

1.00 

78.9070739 

1.3 

6.3 

1.  10 

86.7977809 

-105.9 

-107. 1 

1.20 

94.6884870 

-169.3 

-174.8 

Figure  3.  5.     Pressure  as  a  fxmction  of  volume.     The  solid  curve  is  for 
pressures  determined  using  the  virial  theorem  and  the 
dashed  curve  is  for  pressures  determined  using  the 
derivative  of  the  energy  vs  volume  curve  at  V_.     (Vq  is 
the  volume  corresponding  to  the  experimental  lattice 
parameter. ) 
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first  derivative  of  the  pressure  as  a  function  of  volume  evaluated  at 

V  ,  using  Eq.    (2.47)  as  described  in  Sec.   2.  5.     The  bulk  modulus  was 
0 

also  determined  as  the  volume  times  the  second  derivative  of  the 
total  energy  as  a  function  of  volume  evaluated  at  V  ,   as  given  in 
Eq.    (2.  48).     This  value  of  the  bulk  modulus  was  the  same  as  that 
determined  from  the  slope  of  the  pressure  vs.   volume  curve 
resulting  from  the  first  derivative  of  the  energy  as  a  function  of 
volume.     The  only  differences  in  these  tv/o  va].ues  of  the  bulk 
modulus  were  attributed  to  the  inaccuracies  in  the  numerical  methods 
used  to  evaluate  the  first  and  second  derivatives,   as  would  be 
expected. 

The  compressibility,  which  is  the  reciprocal  of  the  bulk 

-4        -1 
modulus,  was  determined  to  be  6.  96  x  10       Kb       for  the  values  of 

-4        -1 
pressure  obtained  using  the  virial  theorem,   and  7.  16  x  10       Kb 

for  the  values  of  pressure  obtained  using  the  total  energ;>'  as  a  function 

of  volume.     These  two  values  differ  from  the  experimental  value 

of  7.  49  X  lo"     ?:b'    ,  by  7%  in  the  first  case  and  4%  in  the  latter  one, 

which  is  reasonably  good  agreement. 

3.4    Energ)^  Bands  and  Density  of  States 
The  energy  bands  for  five  calculations  of  copper  for  the  lattice 
parameters  corresponding  to  the  experimental  value  and  10%  and  20% 
volume  expansion  and  compression     are  given  in  Figures  3.  6,   3.7, 
3.  8,    3.  9,   and  3.  10.     The  energy  bands  in  these  figures  are  shown 


Figure  3.  6  (a  &  b)    Energy  bands  for  FCC  copper  with  a  =  0.  7225 
and  A  =  6.  321  a.u.     (c)  Density  of  states  for  this 
calculation. 
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Figure      3.6  a 
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Figure        3.6b 
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Figure  3.  7.     (a  &  b)    Energy  bands  for  FCC  copper  with  a  =  0.  7225 
and  A  =  6.  574  a.  u.     (c)  Density  of  states  for  this 
calculation. 
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Figure  3.  8.     (a  &.  b)    Energy  bands  for  FCC  copper  with  Q!=  0.  7225 
and  A  =  6.  809  a.  u.     (c)  Density  of  states  for  this 
calculation. 
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Figure     3.8b 
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Figure  3.  9.     (a  &  b)    Energy'-  bands  for  FCC  copper  with  a  =  0.  7225 
and  A  =  7.028  a.u.     (c)  Density  of  states  for  this 
calculation. 
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Figure      3.9a 
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Figure      3.9  b 
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Figure  3.  10.     (a  &  b)    Energy  bands  for  FCC  copper  with  q:=  0.  7225 
and  A  =  7.  235  a.  u.     (c)  Density  of  states  for  this 
calculation. 
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Figure      3.10a 
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Figure     3.10  b 


67- 


a: 


O 

c 


3 


■68- 


in  the  directions  of  higliest  symmetry  in  the  first  Brilloiiin  zone  of 

the  FCC  structure.     The  value  of  a  used  for  all  five  calculations 

was  0.7225.     The  corresponding  density  of  states  for  each  calculation 

is  also  given. 

Energy  differences  for  these  calculations,    that  represent  the 

bandwidths  and  relative  locations  of  the  bands  with  respect  to  each 

other  and  with  respect  to  the  Fermi  energy,  are  given  in  Table  3.  4. 

The  energy  differences  given  as  (T    ^'  ~r^)  and  (X    -  r    )  give  the 

relative  location  of  the  d-band  witli  respect  to  the  bottom  of  the  sp- 

band,  while  the  differences  (E  -X^)  and  (E„  -  L„)  locate  the  top  of 

ID  to 

the  d-band  with  respect  to  the  Fermi  energy'.     The  energy'  differences 
(Xp  -  X  )  and  (X/  -  r  ')  give  the  bandwidths  of  the  d-band  and  sp- 
band,   respectively.     (E    -  L    ')  represents  the  location  of  the  high 
symmetry  point  with  an  eigenvalue  nearest  the  Fermi  energy.     This 
energy  difference  also  gives  a  measure  of  the  "neck"  radius  of  the 
copper  Fermi  surface. 

A  similar  set  of  energy  differences  for  experimental  results, 
results  of  previoxis  calculations,  and  results  of  present  calculations 
for  different  values  of  o.  (at  the  equilibrium  lattice  constant)  are  given 
in  Table  C.  ]  in  Appendix  C,  and  are  not  included  in  this  table. 

The  relative  movements  and  changes  in  bandwidth  of  the  bands 
as  the  lattice  parameter  is  increased  are  as  expected;  namely,  both 
the  d-band  and  the  sp-band  narroAv  and  the  d-band  moves  toward  the 
bottom  of  the  sp-band.     If  one  were  to  continue  increasing  the  lattice 
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parameter,   the  d-band  would  continue  to  become  narrower  and 
eventually  drop  completely  below  the  sp-band,  approaching  the  3d 
atomic  state  in  the  limit  of  infinite  lattice  parameter.     The  Fermi 
energy  "follows"  the  d-band,   moving  tov/ard  the  bottom  of  the  sp- 
band  faster  than  does  the  energy  of  the  L   '  point,   resulting  in  a 
decrease  in  the  energy  difference  (E  -L    ')  and  a  decrease  in  the 
"neck"  radius  of  the  copper  Fermi  surface. 

In  comparing  the  results  given  in  Table  3.  4  v/ith  experimental 
results,   only  those  for  the  experimental  lattice  parameter  (V/Vq=1.0) 
can  be  considered.     For  these  results,  the  d-band  appears  to  be 
too  broad  and  too  close  to  the  Fermi  energy  when  compared  with  the 
results  of  photoemission  studies.'  The  (E^-L^')  energy  difference 

also  appears  too  large  when  compared  to  the  results  of  Spicer'      or 

39 
Lindau  and  Walldon.         However,   the  value  of  such  comparisons  is 

questionable.     Nevertheless,   the  value  of  (E  -L^'^  given  here  as 

0.  063  Ry  is  only  slightly  larger  than  the  value  of  0.  061  Ry  given 

by  Janak  et  al.        as  being  the  value  that  gives  a  "neck"  radius  in 

agreement  with  experiment. 


CHAPTER  IV 

CONCLUSIONS 

The  main  purpose  of  this  dissertation  was  the  evaluation  of 
the  Xa  exchange  approximation  as  used  in  self- consistent  APW 
calculations  for  a  metal,   such  as  metallic  copper.     To  achieve  this, 
the  total  energy  as  a  function  of  lattice  parameter  was  determined 
and  these  results  were  used  to  determine  cohesive  energy'  and 
compressibility.     Comparison  of  these  quantities  with  experimental 
results  gives  a  measure  of  the  accuracy  of  these  calculations  and 
some  indication  of  the  worth  of  the  methods  employed.     Of  course, 
a  complete  evaluation  of  the  methods  would  involve  a  large  number 

of  calculations  on  a  variety  of  materials;  such  other  calculations 

26    27  28 

have  already  been  done  by  Averill      '        and  Hattox. 

To  evaluate  the  results  of  the  present  work,   we  first  consider 

the  method  proposed  to  determine  the  Q;  that  was  used  in  the  Xo: 

exchange  approximation.    As  was  pointed  out  in  Sec.    3.  1,   the  value 

of  o:  that  woidd  give  zero  pressure  at  the  experimentally  determined 

lattice  spacing  was  found  to  be  0.  7225.     The  results  given  in  Table  3.  3 

show  that  the  calculated  pressure  at  the  experiinental  lattice  spacing 

resulting  from  the  use  of  this    value  of  a  was  1.  3  Kb  instead  of  zero 

as  required.     However,   this  represents  a  difference  in  the  potential 
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energy  and  twice  the  kinetic  energy  of  approximately  0.002  Ry. 
This  is  a  very  small  fraction  (on  the  order  of  3  x  10      )  of  the  6562  Ry, 
which  is  the  approximate  magnitude  of  these  two  nuinbers.  The  virial 
coefficient,   defined  in  Eq.    (2.  20)  as  the  ratio  of  the  negative  of  the 
potential  to  twice  the  kinetic  energy,   was  0.  9999997  for  this  calculation, 
or  very  nearly  1.  0,   the  desired  value. 

The  total  energy  as  a  function  of  lattice  parameter,   as  one 
would  expect  for  this  calculation,    has  a  minimum  very  close  to  the 
experimentally  determined  lattice  parameter,   as  can  be  seen  in 
Figs.    3.  3  and  3.  4.     The  pressure  of  6.  3  Kb, given  in  Table  3.  3, 
indicates  that  the  minimum  occurs  slightly  outside  the  experimental 

value.     This  is  evident  if  one  considers  that  6.  3  Kb  represents  a 

-  5 
slope  in  the  energy  vs.   volume  curve  of  about  -4x10       in  atomic  units. 

The  cohesive  energy  obtained  from  the  calculation  at  the  experi- 
mental lattice  spacing  with  0:=  0.  7225  was  0.  286  Ry.     This  value  is 
within  11%  of  the  experimental  value  of  0.  2  57  Ry,   and  is  probably 
within  the  experimental  accuracy  of  determining  tlie  cohesive  energy 
at  O^K. 

The  energy  bands  for  this  calculation  agree  reasonably  well 
with  the  results  of  photoemission  studies  given  in  Table  C.  1,    though 
such  direct  comparisons  are  somewhat  questionable.     However,   the 

resulting  (E  -L    ')  energy  difference  of  0.  063  Ry  agrees  very  well  with 

40 
the  0.  061  Ry  given  by  Janak  et  al.        as  being  the  \^lue  that  corresponds 
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to  a  "neck"  radius  of  the  copper  Fermi  surface  in  agreement  with 
experiment. 

Probably  the  most  significant  results  achieved  in  the  present 
calculations  were  the  agreement  between  the  two  sets  of  calculated 
values  of  the  pressure  as  a  function  of  lattice  parameter  and  the  values 
of  the  compressibility  determined  from  these  tv/o  sets  of  data.     The 
agreement  between  these  two  sets  of  calculated  pressures,   given  in 
Table  3,  3,   is  evident  in  Fig.   3.  5.     The  compressibilities  deterinined 
from  the  slope  of  these  two  curves  agree  with  the  experimental  value 
(to  within  7%  and  4%). 

Hence,   one  can  conclude  that  the  XOi  method,  as  applied  in  these 
calculations,   can  be  used  to  obtain  total  energy  as  a  function  of  lattice 
parameter,   cohesive  energy,   and  com.pressibility  that  are  in  good 
agreement  v/ith  experimental  values.     However,   these  results  would 
have  agreed  with  experiment  just  about  as  well,   if  the  a.  that  satisfied 
the  virial  theorem  for  tlie  atomic  calculation  (a=  0.  70635)  had  been 
used  in  the  present  calculations,   even  though  the  minimum  of  the 
energy  vs.  lattice  pai-ameter  curve  wovild  be  slightly  outside  the 
experimentally  determined  lattice  parameter,  by  approximately  0.  5%. 
This  indicates  that,   from  a  practical  point  of  view,   this  value  of  a.  is 
probably  the  best  one  to  use  in  a  calculation  of  this  type,   as  was  the 
case  for  the  calciilations  by  Averill     *        and  Hattox. 


APPENDICES 


APPENDIX  A 


METHOiOS   OF  IMPROVING  THE  ACCURACY 

OF  THE  LOGARITHMIC  DERIVATIVES 

USED  IN  AN  APW  CALCULATION 


A.  1     Introduction 

In  using  the  APW  method  for  calculating  the  energy-bands  of  a 
crystal,  the  determination  of  pressure  and  cohesive  energy  requires 
greater  precision  than  is  needed  for  determining  just  the  energy 
bands  E(k).     There  are  several  obvious  v/ays  of  improving  accuracy, 
such  as  increasing  the  number  of  augmented  plane  waves  in  the 
expansion  of  the  wave  functions  or  increasing  the  number  of  spherical 
harmonics  used  to  expand  each  a.ugmented  plane  wave,   or  both.     The 
potential  function  itself  can  be  improved  to  include  "non-muffin  tin" 
terms. 

A  more  fundamental  method  of  improving  accuracy  is  to  improve 
the  numerical  methods  used  to  solve  the  radial  Schrodinger  equation, 
given  as 

for  the  radial  function  U     and  the  numerical    methods  used  to  evaluate 
U    (Rs)/U    (Rs),   the  logaritlimic  derivatives  that  appear  in  the  matrix 
elements  of  the  secular  equation.     The  improvements  to  the  methods 
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used  in  the  present  work  will  be  discussed  in  detail  in  the  following 
sections. 

A.  2    Starting  Values 
To  solve  the  radial  Schrodinger  equation,   Eq.   (A.  1),   one  solves 


for  P    (r)  given  by 


P    (r)  =  r  U    (r).  (A.  2) 


The  differential  equation  then  becomes 

p'^'(r)  =  g(r)  P^(r)  (A.  3) 

where  P  "(r)  represents  the  second  derivative  of  P  .  (r)  v/ith  respect 
to  r,  and  g(r)  is  given  by 

g(r)  =  -E-  V(r)+-^^^        .  (A.  4) 

To  solve  Eq.    (A.  3)  numerically,   a  mesh  of  radial  points  is  set  up 

with  the  value  of  the  potential  V(r  )  given  for  each  point  n  of  the  mesh. 

The  selection  of  this  mesh  is  discussed  in  detail  in  Sec.  A.  4.     The 

Noumerov^'^  method  is  then  used  to  solve  for  P /r)  for  each  value  of  £ 

to  bo  used.     This  method  is  used  to  dctermjine  the  value  of  P^  (r^) 

from  the  value  of  E,  V(r  ),   P  (r      J,   and  P    (r        ).     Therefore,  to 

n  j2,   n- 1  X,     ]i~  ^ 

start  this  "outward"  integration  scheme  of  Eq.    (A.  3),  the  starting 
values  of  the  radial  function,   P/j^J  and  IP    (r^),   must  be  supplied. 
Since  the  accuracy  of  the  final  wave  function  is  dependent  on  the 
accuracy  of  approximating  these  two  starting  values,   they  must  be 
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determined  with  some  care.     Greisen     has  shown  that  the  approximation 

P    (r  )^   r    ^"^^  (A.  5) 

in—      n 

can  result  in  errors  in  the  eigenvalues  of  the  APVv''  calculation  that  are 

large  compared  to  the  self- consistency  criteria;  these  errors  are 

expecially  large  for  functions  with  a  large  amount  of  s-like  symmetry. 

This  section  is  devoted  to  the  development  of  a  more  elaborate  scheme 

of  approximating  these  starting  values,   similar  to  those  given  by 

Hartree. 

First  the  potential  is  expanded  in  a  power  series  of  r 

00 

r   -V     =    y^     b     r   "^  (A.  6) 

n      n       ^— <^      m  n 

m-0 

where 

V      =    V(r  ). 
n  n 

For  small  r,   the  expansion  need  be  carried  to  only  a  few  terms, 

such  as 

2 

r   'V    ~b     +b   r     +b   r  ,.       , 

n      n        0         In        2  n  (A.  7) 


v/here 


^im^(r   .V   )    =   -  2Z    =  b^  (A.  8) 

r  ^0     n      n  0 


and  where  Z  is  the  nuclear  charge.     To  solve  for  the  other  tvv^o 
coefficients,    Eq,    (A.  7)  is  evaluated  for  n  =  1,    2,   yielding 

r^.V^^    -2Z  +  b^r^  ^-b^r^^  (A.  9) 


and 


rg-y^-    -2Z  +  b^r2  +b2r22      .  (A.  10) 


-78- 


Equations  (A.  9)  and  (A.  10)  can  then  be  solved  for  b^  and  b^,  yielding 

(r  V     -  r   V  )  2Z(r     +r  ) 

b     -_-^ -i-^     +    i ^  (A.  11) 

1        (r2  -  r^)  r^   r^ 


and 


(V   -V   )  „ 

b     =7-^— T-      -     — -^  .  (A.  12) 

2      (r^-r^)  r^r^ 

Then  P    (r  )  is  expanded  in  a  power  scries  given  by 
in 


^+1  Y, 


m 


P    (r  )    =    r    ^"'      ^.    a     r  "'  (A.  13) 

I     vi  n  m=0      m  n 

v/here  a    is  an  arbitrary  constant.     To  approximate  Pjj(rj)  and  V ^{v^, 

only  the  first  four  terms  of  the  expansion  are  used,  yielding 

P    (r  )  ~  a   r  +  a,r  +  a„r  +  a„r  .  (A.  14) 

X^n-On  In  2n  3n 

Substituting  the  expansions  of  ^  J<^^  and  r^.  V^,   given  in  Eqs.    (A.  14) 

and  (A.  7),   into  the  differential  equation,   Eq.    (A.  3),   the  a's  can  be 

determined  by  requiring  that  the  coefficients  of  the  power  of  r  to 

vanish  independently  yielding 

a     ---^Q^O—  (A.  15) 

^1  "       2(X  +  1) 

a.(E+b J  +  a  b 
_  __Q 1 1_Q.  (A.  16) 

2  2(2j£+3) 


and 


s^ -^  V2 -^  ^i<^"N^  ^'^'^'^ 


^3       '  6  (X+  2) 


With  these  values  of  a^,   2.^,   and  a^,   and  setting  (the  so-far  arbitrary) 
a     such  that  the  value  of  the  resulting  radial  function  at  any  point  never 
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exceeds  the  numerical  limits  of  the  computer,   the  resulting  approximation 
to  the  starting  values  gives  very  accurate  eigenvalues  for  the  APW 
calculation. 

The  Noumerov  method  itself  can  be  improved  by  a  method  given 

44 
by  Fischer,  but  it  v/as  felt  that  this  correction  vrould  be  small 

compared  to  that  resulting  from  the  improved  starting  values  and  has 

not  been  employed  in  the  present  calculation. 


A.  3    Evaluation  of  the  Derivatives  at  R 


The  logarithmic  derivatives,   P  ,   given  by 


U    (r) 

D      = — 

l       U^(r) 


(A.  18) 
r=R 


s 

(where  R     is  the  radius  of  the  APW  sphere)  can  be  evaluated  using 
s 

I 

P  '(r)  l^^pj     =  P      (R  )  by  the  relationship 

P  '  (R  )  P   (R   ) 

U    '  (R  )  =       ^^      ^       -  -^-2 (A.  19) 

is  R  ^  R 

s  s 

The  derivative  of  P    (R  )  can  be  evaluated  very  accurately  if  the  fact 

Jo  S 

is  used  that  P  (r)  is  a  solution  to  the  differential  equation 
P^'(r)  =  g(r)  P^(r) 

given  as  Eq.    (A.  3)  above. 

First,   the  equations  for  P(n+1)  and  P(n-l)  are  given  as 

2  3 

P(n+1)  =  P(n)  +  hp'(n)  +~  p"(n)+^  p"'(n)  +  ...  (A.  20) 
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and 

P(n-l)  =  P(n)  -  hp'(n)+f7    p' '(n)  -  rp  P'    (n)  +  . . .  (A.21) 


h^        '\   X      iL        '" 


21  "   '      3! 

where 

Pfa)  .  P/r_^) 

and  h  is  the  spacing  of  the  radial  mesh  points.       Subtracting  (A.  21) 

from  A.  20  yields 

3  5 

P(n+1)  -  P(n-l)  =  2hp'(n)  +J^  p'"(n)  +~  P^(n)  +  ....  (A.  22) 

Similar  equation  for  P(n+2)  and  P(n-2)   are  then  set  up,  the  difference 
of  which  is  given  by 

P(n+2)  -  P(n-2)  =  4hp'(n)  +^  p"'(n)  -^^f^  P^(n)  +  . . . .  (A.  23) 

Multiplying  Eq.    (A.  22)  by  8  and  subtractmg  Eq.   (A.  23)  from  it  yields 


8  rP(n+l)  -  P(n-l)]    -    [p(n+2)  -  P(n-2)J 


=  12hp'(n)-^p''(n)  +  ....  (A.  24) 

Similarly  Eq.    (A.  22)  is  multiplied  by  2  and  subtracted  from  Eq.    (A.  23), 
yielding 

rP(n+2)  -  P(n-2)1     -2  rP(n+l)  -  P(n- l)j 

,     l|k%'"(,^),-60|!    pV(,)^....  (A.25) 

Taking  the  second  derivation  of  Eq.    (A.  25)  and  using  the  differential 
equation,   Eq.    (A.  3),   gives 

rg(n+2)P(n+2)  -  g(n-2)P(n-2)]    -     2  f  g(n+l)P(n+l)    -     g(n-l)P(n-l)J 

.  JJ^    p  V         ^  SmL  j.^%)  ^  .  . .  .  (A.  26) 

?,  I  ^  5! 
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This  equation  is  then  multiplied  by  "  ,   solved  for      r^  P    (n)  and 
substituted  into  Eq.    (A.  24).     Dropping  the  terms  involving  the  seventh 
and  higher  derivatives,  and  solving  for  P'(n)  yields  the  final  relationship 

given  by 

12  2 

[l  -—  g(n-2)]p(n-2)  -  [s  -  ^  g(n- 1)]  P(n- 1) 

2  2  1 

+  [s  -  ^  g(n+l)]  P(n+1)    -  [l-  ~  g(n+2)]   P(n+2)  (A.  27) 

which  was  used  in  Eq.    (A.  19)  to  evaluate  the  logarithmic  derivatives 
in  these  calculations. 

A.  4    Effects  of  the  Choice  of  Radial  Mesh 
In  the  application  of  numerical  methods  to  the  solution  of  the 
radial  Schrodinger  equation,   Eq.    (A.  3),   a  radial  mesh  must  be  set  vi.p 
for  a  finite  number  of  points,  usually  between  100  and  1000  points. 
The  spacing  of  the  radial  points  must  be  small  at  small  values  of  the 
radius  (r)     to  insure  accurate  results  in  the  region  v/here  both  the 
potential  function  and  the  wave  function  are  varying  rapidly.     The 
spacing  must  then  be  increased  (if  we  are  to  keep  the  total  number 
of  points  within  reason)  as  r  increases  so  that  the  mesh  will  cover 
the  range  of  r  required  for  the  calculation. 

There  are  essentially  two  ways  of  defining  such  a  radial  mesh. 
The  first,   known  as  the  logarithmic  mesh,   is  given  by 

r     =  R^e""^  (A.  28) 

n         0 

where  R     and  h  are  set  to  satisfy  the  spacing  and  range  requirements 
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of  the  mesh,   as  mentioned  above.     The  main  advantage  of  this  type 
mesh  is  its  simplicity.   However,   as  will  be  sho^^TJ  later,    this  is  also 
the  main  disadvantage.     The  second,   called  a  linear  mesh,   is  set  up 
in  blocks  of  mesh  points,   in  which  the  mesh  spacing  is  constant 
betv/een  each  mesh  point.     In  this  type  of  mesh,  the  spacing  is  set  up 
to  be  small  in  the  first  block  and  is  doubled  as  r  goes  into  each 
succeeding  block.     These  blocks,   usually  less  than  10  in  number,  need 
not  contain  the  same  number  of  points  and  can  be  adjusted  so  that  the 
range  of  r  values  required  is  covered  while  retaining  some  control  of 
the  spacing  at  the  extremes  of  r.     This  versatility  is  the  main  advantage 
of  the  linear  mesh  and  the  lack  of  it,   the  main  disadvantage  of  the 
logarithmic  mesh. 

To  test  the  accuracy  that  can  be  achieved  when  solving  the  radial 
Schrodinger  equation  numerically  for  these  two  types  of  radial  meshes, 
the  logarithmic  derivatives  evaluated  at  R    were  determined  for  a  pure 
coulombic  potential  -2Z/r  by  both  numerical  and  analytical  methods. 
These  results  are  tabulated  and  the  compared  for  the  two  radial  meshes 
considered.      Tests  on  the  accuracy  of  evaluating  the  eigenvalues  of  the 
core  states  on  a  linear  mesh  are  also  presented. 

The  logarithinic  derivatives  for  16  different  logarithmic  radial 

meshes  are  given  in  Table  (A.  1)  for  jii=  0,    1,   and  2.     The  energy  was 

set  at  0.  1  Ry,   which  is  representative  of  the  energy-band  states.     The 

values  of  Z=26  and  R  =2.  5  a.u.  were  selected  so    that  a  comparison 
s  ^ 
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TABLE  A,  1 


The  logarithmic  derivatives  U^/U^ ,  evaluated  at  R^  =  2.  5,   for  a  -^Z/r 
potential,  with  Z  =  26,  on  a  logarithmic  radial  mesh.     The  energy  in 
all  cases  was  0.  1  Ry.     Ri  is  the  minimum  r  value  used  and  N  is  the 
total  number  of  radial  mesh  points. 


(a.  u. ) 

N 

AR 

@  R 
s 

^:/". 

(jJ-O)  X  10-2 

(2  =  1)  X  10-2 

(i  =  2)   X  10-1 

10-2 

200 

0.071 

4. 59011 

1.  15084 

3.96734 

10-2 

300 

0.047 

4. 59583 

1. 15117 

3.  96782 

10-=^ 

400 

0.035 

4. 59684 

1.  15123 

3.96790 

10"2 

500 

0.028 

4. 59716 

1. 15124 

3.96792 

10-* 

200 

0.  132 

4. 53063 

1. 14639 

3.  96066 

10-* 

300 

0.087 

4. 59258 

1.  15033 

3.96662 

10-* 

400 

0.065 

4. 60286 

1.  15097 

3.  96753 

10-* 

500 

0.051 

4.60  563- 

1.  15114* 

3. 96777* 

10-^ 

200 

0.  194 

4.27631 

1.  12869 

3.93144 

io-« 

300 

0.  127 

4.  54079 

1.  14705 

3.96168 

10-^ 

400 

0.09  5 

4.  58566 

1.  14996 

3.  96607 

10-^ 

500 

0.075 

4. 59900 

.1.  15073 

3.  96719 

10-  = 

200 

0.258 

3.72558 

1.08247 

3.84538 

10-^ 

300 

0.  168 

4.41285 

1. 13847 

3.94797 

10-^ 

400 

0.  125 

4.  54544 

1.  14735 

3.96214 

10-^ 

500 

0.099 

4. 58218 

A . — • 

1.  14968 

3.  96566 

*Value  nearest  analytic  results. 
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co'Lild  be  made  with  analytical  results  made  available  by  Wood. 
The  results  of  this  table  are  displayed  graphically  in  Figs.   (A.  1), 

(A.  2),   and  (A.  3).     The  main  feature  of  interest,   which  shows  up 

-4 
especially  well  in  Fig.    (A.  1),   is  tha.t  R     ^^  10       a.  u. ,   the  accuracy 

of  the  logarithmic  derivatives  is  dependent  almost  entirely  on  the  mesh 

spacing  at  R   .     This  is  also  true  in  Figs.    (A.  2)  and  (A.  3).   but  it  is 

not  as  obvious.     This  does  show  up  in  Table  (A.  1)  in  all  three  cases, 

where  the  most  accurate  values  of  the  logarithmic  derivatives  are 

-4 
achieved  with  500  points  and  R     =  10      .     The  relative  errors  for 

these  three  results  are  0.0004,   0.0001,   and  0.00004  for  l^  0,    1,   and 

2,  respectively. 

The  relative  err'ors  in  the  logarithmic  derivatives,   evaluated 

for  X=  0,    1,   and  2,   on  four  linear  meshes  are  given  in  Table  (A.  2), 

and  the  relative  errors  in  the  eigenvalues  of  the  "core"  states  are 

given  in  Table  (A.  3)  for  the  same  four  linear  meshes.     It  should  be 

-4 
pointed  out  tliat  not  only  is  R  ~  10       in  all  four  cases,  but  the  spacing 

is  so  constructed  as  to  maintain  a  smaller  spacing  at  R     than  can  be 

achieved  with  the  logarithmic  mesh. 

Comparing  the  results  for  a  200  point  linear  mesh  with  the 

results  for  the  "best"  500  point  logarithmic  mesh,   the  relative  error 

is  one  tenth  as  large  for  the  linear  mesh.     It  can  easily  be  understood 

why  the  linear  mesh  was  selected  for  use  in  the  present  calculations. 
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TABLE  A.  2 


The  relative  error  in  the  logarithmic  derivatives  for  four  linear  radial 
meshes,  v/ith  Rg  =  2.  5  a.  u. ,   Z  =  26  a.  u. ,  and  E  -  0.  10  Ry,   for  a 
-2Z/r  test  potential.     Rj  is  the  minimum  value  of  r  used  and  N  is  the 
total  number  of  points  in  the  radial  mesh. 


N 

R    X  104 
\'a.u.) 

AR 
©Rg 

Relative  error  in  U'  /U 

;.=  0       1       X-  1 

1=  2 

200 

2.6 

0.017 

6.5  X  lO"^ 

1.4  X  lO"^ 

1.  1  X  10"5 

300 

2.  1 

0,013 

8.2  X  10-6 

3.  5  X  lO"^ 

2.1  X  lO"^ 

400 

1.7 

0.011 

5.  3  X  10"^ 

1.  5  X  10"^ 

7.0  X  lO"'^ 

500 

1.4 

0.009 

2.  9  X  10-6 

7.0  X  lO"'^ 

3.0  X  lO"*^ 

i 
J 
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TABLE  A.  3 


Relative  error  in  the  core  state  eigenvakies  for  four  linear  meshes, 
with  Rg  =  2.  5  a.u. ,   Z  =  26  a.  u. ,  and  a  -2Z/r  test  potential.     R]^  is 
the  minimiun  value  of  r  used  and  N  is  the  total  number  of  points  in 
the  radial  mesh. 


' 

N-200 
R  =2.6  X  10-4 

N  =  300 
Rj  =  2.  1x10'^ 

N  =  400 
R   =1.7x  lO' 

N  =  500 
R  =1.4  X  lO" 

ARj;.    =  0.017 
s 

ARr    =0.013 
^s 

ARj^   =  0.011 

s 

ARr    =  0.009 
^s 

— 

Is 

2.4  X  lO"^ 

3.9  X  lO''^ 

9.  9  x  lO" 

2.  1  X  lO"^ 

2s 

1.8  X  lO"^ 

1.7  X  lO''^ 

2.9  X  lO'^ 

1.3  X  lO"^ 

2p 

-6 
1.0  X  10 

1.  1  X  lO''^ 

2.2  X  lO"^ 

1.0  X  lO"^ 

3s=:= 

1.8  X  lO" 

-6 
4.  1  X  10 

1.8  X  lO' 

6.9  X  lo""^ 

3p^= 

1. 5  X  lO' 

2.6  X  lO'^ 

1.2  X  lO" 

5.  2  X  lo'"^ 

^These  states  are  considered  as  band  states  in  the  present  calciilations. 


APPENDIX  B 


TOTAL  ENERGY  CALCULATIONS  USED  TO  EVALUATE  THE 
"FROZEN  CORE"  MODEL 


The  "frozen  core"  model  of  a  self-consistent  APW  calculation  is 

based  on  the  assumption  that  the  core  state  functions  do  not  change 

much  while  going  from  the  free  atom  to  the  self- consistent  crystal 

result,   and  can,   therefore,  be  held  fixed  throughout  the  calculation 

without  significantly  affecting  the  results.     Some  insight  into  the 

validity'  of  this  assumption  can  be  obtained  if  the  total  energy  is  broken 

up  such  that 

E  =  E         +  E         +  E  (B.  1) 

c-c         v-v         c-v 

where  E         is  that  part  of  the  total  energy  resulting  from  the  kinetic 

energy  of  the  core  electrons  and  that  part  of  the  potential  energy 

i^epresenting  the  coulomb  and  exchange  interactions  of  the  core  electrons 

with  the  nucleus  and  the  other  core  electrons,     E         represents  the  same 

v-v 

energy  for  the  valence  electrons,   and  E         represents  the  interactions 

between  the  core  and  valence  electrons.     The  E         energies  for  the 

c-c 

self- consistent  crystal  calculations  and  the  free  atom  calculation  can 
be  determined  and  compared,   giving  an  indication  as  to  hov/  these  core 
states  vary,   and  as  to  what  effect  these  variations  might  have  on  total 
energy  calculations. 
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To  evaluate  these  energies,   the  total  charge  density  is  divided 
into  that  associated  with  the  core  states,    P^M.   and  that  associated 
with  the  valence  state,    p    (r),   such  that 

p{r)  =   p^(r)+  p^(r).  (B.  2) 

The  total  core-core  energy,   E^_^,   is  then  given  by 

where 

r  r  rl^c    P    ^1^') 

-^+-^  o    fr'kir-  -  2        ^-^—  dr'+C  (B.4) 


V^>  =  "    r     "r  Pe(^'>^-'-2,  r- 


-+■;    rp^(r')dr--2j 

and  Sv  ^  is  the  total  kinetic  energy  of  the  core  electrons.     The  total 
c    i 

valence-valence  energy,   E   _    ,   is  given  by  the  same  four  equations 
with  the  subscript  (c)  replaced  by  the  subscript  (v).     The  energy 
associated  with  the  core- valence  interactions  is  given  by 


0  r 

h/P    (r)\    i  (B.5) 


V.     =-  2  /  ^^s-^^  dr'  +  C  ■  (B.6) 


E  =  E  -  (E         +E        ).  (B.7) 

c-v  c-c         v-v 


The  results  for  a  calculation  on  copper  with  Oi=  0.7225,   and  two 
core  configurations  are  given  in  Table  (B.  1).     The  most  interesting 
results  are  the  differences  in  core- core  energies  for  the  tv/o  configura- 
tions.    If  the  configuration  of  the  "frozen  core"  was  selected  to  be 
ls^2s^2p^,    (ls-2p),   the  difference  of  the  core-core  energy  would 
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TABLE  B.  1 


Total  enex'gies  for  self- consistent  atomic  and  crystal  calculations  on 
copper  with  a=  0.7225,   for  two  core  configurations. 


1 
Energy  (Ry) 

Crystal 

Atomic 

Difference 

(ls-3p)  Core 
Core-Core 
Val.  -Val. 
Core-Val, 

-3177.7552 

-    552.4849 

449.  1419 

-3177. 5900 

-    550. 5609 

447.3383 

0.1653 

1.9240 

-1.8036 

(ls-2p)  Core 
Core-Core 
Val.  -Val. 
Core-Val. 

-2885.6566 

-    938.8G40 

543.4224 

-2885.  6624 

-    937.3516 

542.2015 

-0.0058 

1. 5124 

-1.2209 

Total 

-3281.0982 

-3280.8125 

0.28  57* 

=''Cohesive  energy. 
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indicate  that  this  should  be  an  acceptable  model.     However,   if  the 

o         o         f^         O         fj 

configuration  is  chosen  to  be  Is   2s   2p   3s   3p  ,    (ls-3p),  which 
is  usually  done,   the  difference  in  core-coi-e  energy  is  on  the  order 
of  the  cohesive  energy  which  indicates  changes  in  some  of  these 
core  states,  namely  the  3s  and  3p  states,   that  could  have  significant 
effects  on  such  quantities  as  cohesive  energy  or  pressure,   and  the 
use  of  such  a  "frozen  core"  model  would  be  somev/hat  limited. 


APPENDIX  C 


THE  ENERGY  BANDS,    DENSITY  OF  STATES, 

TOTAL  ENERGY,  AND  COHESIVE  ENERGY 

OF  COPPER  AS  A  FUNCTION  OF  a 


C.  1    The  Energy  Bands  and  Density  of  States 
This  section  is  devoted  to  the  presentation  of  the  energy  bands 
and  density  of  states  for  the  four  values  of  a.  used  to  determine  the 
a.  that  yields  zero  pressure  at  the  observed  lattice  spacing,   as  described 
in  Sec.    3.  1.     Each  calculation  was  done  for  2048  points  in  the  first 
Brillouin  zone  and  used  a  500  point  linear  radial  mesh  (for  integration 
of  the  radial  Schrbdinger  equation).      Though  the  self- consistency 
criterion  used  in  these  calculations  was  applied  to  the  pressure  variation 

between  two  successive  iterations,   the  variation  of  any  eigenvalue  between 

-4 
tv/o  successive  iterations  was  found  to  be  less  than  10       Ry  for  the  core 

-  5 

states  (usually  the  Is  state  has  the  largest  variation)  and  less  than  10 

Ry  for  the  band  states.     The  resulting  change  in  total  energy  betv/een 
two  successive  iterations  was  about  10       Ry. 

Figures  C.  1,   C.  2,   C.  3,   and  C.4  show  the  energy  bands  for 
copper  for  four  values  of  a,   along  the  directions  of  highest  symmetry 
in  the  first  Brillouin  zone.     All  four  correspond  to  an  APW  sphere 

radius  (R  .)  of  2.  407  a.  u. ,  which  is  one  half  the  nearest  neighbor 

s 

distance  for  the  experimentally  determined  lattice  spacing  extrapolated 


■97- 


Figure  C.  1.     (a&b)    Energy  bands  for  FCC  copper  with  a  =  0.  700  and 

R  ,  ^  2.407  a.\i.     (c)  Density  of  states  for  this  calculation. 
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Figure     C . 1  a 
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Figure      C.l  b 


101- 


LLi     U_ 


(3)N 


Figure  C.  2.     (a&b)    Energy  bands  for  FCC  copper  with  a  =  0.  70635 
and  R     =  2.407  a.  u.     (c)  Density  of  states  for  this 
calculation. 
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Figure      C.2a 
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Figure      C.2b 
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Figure  C.  3.     (a&b)    Bnergsr  bands  for  FCC  copper  with  a  =  0.  720  and 

R     -  2.407  a.  u.     (c)  Density  of  states  for  this  calculation. 
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Figure       C.3b 
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Figure  C.4.     (a&b)    Energ>'  bands  for  FCC  copper  with  0:==  0.  770  and 

Rg  =  2.  407  a.  u.     (c)  Density  of  states  for  this  calculation. 
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Figure     C.4q 


112- 


Figure      C.4b 
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to  0°K.     The  corresponding  density  of  states  for  each  calculation 
is  also  given. 

Table  C.  1  gives  energy  differences  for  these  calciolations, 
that  represent  the  bandwidths  and  relative  locations  of  the  bands  with 
respect  to  each  other  and  with  respect  to  the  Fermi  energy.     The 
energy  differences  (r„j.,  -  r  J  and  (X^-  r  J  give  the  relative  location 
of  the  d-band  with  respect  to  the  sp-band,   while  the  differences 
(E    -  Xj.)  andCE    -  L   )  locate  the  top  of  the  d-band  with  respect  to  the 
Fermi  energy.     The  energy  differences  (X     -  X  )  and  (X      "   T  •, ) 
give  the  bandwidths  of  the  d-band  and  sp-band  respectively.    (E  -L     ) 
represents  the  location  of  the  high  symmetry  point  \\n.th  an  eigenvalue 
nearest  t]ie  Fermi  energy.     This  energy  difference  gives  an  indication 
of  the  "neck"  radius  of  the  copper  Fermi  surface;  the  necks  lie  along 
the  [  111]  directions. 

Also  included  in  the  table  are  energy  differences  resulting 
from  the  calculation  for  a  =  0.  7225,   using  the  experimentally  determined 

lattice  parameter,   and  a  similar  calculation  for  a~  2/3.     (This  latter 

_  3 

calculation  was  carried  to  a  convergence  of  only  10       Ry  for  the 

total  energy. )    Results  of  previous  calculations  and  experiment  are 
also  included  for  comparison. 

The  variation  of  the  relative  Iqcations  and  widths  of  the  bands 
as  a  function  of  a  ^re  as  found  in  previous  calculations.     That  is,   as 
a  is  increased,   the  d-band  narrows  and  its  average  energy  moves 
closer  to  the  bottom  of  the  sp-band  at  V    ,  while  the  sp-band  width 


-115- 


m 


TS 

C 

a 

xn 

ti 

o 

•1-1 

+j 

m 

O 

a 

0) 

+j 

a 

o 

T5 

d 

•  i-i 

. 

-t-> 

■^ 

(Tl 

lO 

^ 

K 

-(-> 

OT 

Ti 

O 

C! 

-^^ 

rt 

^ 

rn 

T1 

;, 

o 

X! 

tM 

C 

tn 

a 

(U 

TS 

o 

fl 

C 

n! 

0) 

^ 

^ 

1 

.<u   a 


b 

bXi 

n 

S-' 

0) 

T5 

n 

•r-l 

W 

!S 

;>^ 


lO 

■>;}<    lO 

CD    CD    05 

lO    CD 

03  CO  ^  r-H  in 

LO 

CO    ^ 

LO    CO    i-l 

O    (M 

CO  c^  t-  CO  -vH 

o 

o  o 

O    O    O 

O    O 

o  o  o  o  o 

o  o  o  o  o 


CD    lO 


(M    CD 
CO    lO 


00  CO 

CO    •<;}< 


r-  ■*  '^ 

O    O    Oj 

CO  00  i~ 


CO    03    •<^< 
O    CD    O 


O    O    O         r- 


O  CD  CO 
O  -^  LO 
CO    CM    CNJ 


cj 

CO 

fl 

(U 

en 

•r-i 

w 

s;' 

O) 

CO 

Oi 

CSl 

to 

CSl 

CO 

CO 

u 

lO 

CS] 

^ 

t— 1 

1-H 

Cvl 

CO 

Lf5 

i 

1—1 

CSl 

1-i 

T-{ 

t-l 

t-l 

T-l 

tH 

o 

o 

o 

o 

o 

o 

o 

o 

(U 

r, 

CO 

M^ 

CO 

CO 

CO 

LO 

^ 

o 

CD 

O) 

CD 

CD 

c:i 

Oi 

CD 

-r-J 

t- 

c~ 

t- 

l>- 

i> 

c- 

t- 

rt 

CO 

co 

ci 

CD 

C3 

o 

CD 

'ps 

o 
I— 1 

u 

lO 

CD 

^ 

o 

CO 

C<1 

CSl 

CD 

rn 

^ 

CO 

CNl 

CO 

irj 

in 

^ 

CSl 

J3 

cs] 

r-i 

CM 

CSl 

CSl 

CSl 

CSl 

CV] 

O 

, 

. 

, 

. 

. 

. 

. 

. 

•.-1 

o 

o 

o 

o 

o 

o 

o 

o 

(U 

?H 

U^ 

en 

CD 

t— 

CO 

o 

o 

lO 

o 

■"vfl     LO     LO 


CO     IT- 
CO   -^ 


tH    CD    CO 
CO    CD    CO 

CO  CO  rfi 


CO    lO    CS!    CD 
lO    LO    LO    LO    ^ 


o  o  o  o  o 


■^  D-  CO  t-  LO 

LO  Tf  CO  T-l  CO 

■sfl  -^  -"^  %f  CO 

O  O  CD  CD  O 


r-l 
I— 1 

fl  o 

O     jU  og 

o    S    S  rt 


cd 


^ 


!> 
a* 

xs 


bJD 


ttn 


0) 

r25     .r-l       CD 


5 


0 


M    ^     rt     g    >    >    > 

t^  W  ^  h 


a    cj 


CO    r-l    O    CO 

o  CSl  r~  00 
c-  i>  t-  CO 


0000 


s  a  £j  a 


-116- 


d 

d 
o 
U 


U 

W 

pq 

<i 


CM 

1 

O    O    CO    CO    CO    C3 
CO    L—    CO    CD    CD    "^ 

o  o  o  o  o  o 

CD  d  CD  d  d  CD 

CO 

1 

CO  (M  '^jf  t-  t-  c:^ 

i-H    CV!    (M    C^J    CM    CO 
T-(    »-l     \rH    T-H    »-l     r-l 

d  d  d  d  d  d 

in 

X 

cs]  crj  o  CO  ^  tn 

O    O     '-I    T-l     rH    CM 

t-<      T-(       T-H       t-(      t— 1      T-l 

<~>    S    (^    CD    Ci    C^ 

1 

X 

CM       »-<       T-1      T-<      T-H      1-1 

o  o  o  o  o  o 

CO    CD    CO    CO    CO    CO 

c>  d  o  ci  o  ci 

J-l 

X 

1 

LO 

X! 

rj<    l^    CO    CO    C<]    CM 

CO  LO  in  in  m  T^i 

CM    CM    CM    CM    CM    C<1 

d  o  o  ci  ci  o 

1 
X 

CO    T-i    CO    1-1    G:)    CO 

t^  CO  in  in  ■>;)<  CM 
in  in  in  in  in  LO 

C)    O    d    Ci    S    CD 

1 
"in 

o  in  CM  CD  in  CO 
CD  ^  ■<d<  o-D  CO   T-t 

^;t^    "^    -*    ^    ^    '^ 

d  d  d  d  d  d 

a=   0.66667 
a=   0.70000 
a=  0.70635 
Ci=   0.72000 
cv=  0.72250 
o;=  0.  77000 

CD    t-    CO    O)    CD    t-    CO    o:) 

cocococo-'^t^-^'^-* 


o        o 

in  CD  'ct' 

<D     (D     <U 


«+H    ^+-^    ^    t*-i    ^    ^+-^    ^-1    ^-< 
<UOCUCD0CUC)0 


CJ  o 

C  d 

0)     O  (U 

Sh     ^  Sm 

0    cu  0 

C+-(      C4_(  C+_( 

O     0)  0) 

p:?  P^  (i; 


117- 


remains  relatively  uneffected.     The  Fermi  energy  follows  the 
movement  and  the  d-band,   moving  toward  r^  as  a  increases. 

The  bandwidths  and  relative  locations  for  the  present  calcula- 
tions are  in  good  agreement  with  those  determined  recently  by  Janak  et  al. 
xistng  the  KKR  method.     The  only  sigTiificant  differences  are  the 
locations  of  the  X^.  and  L„'  eigenvalues  with  respect  to  the  Fermi 
energy.     However,   these  discrepancies  can  probably  be  attributed 
to  the  accuracies  of  the  different  approaches  to  locating  the  Fermi 
energy  in  the  two  different  methods.     Part  of  these  discrepancies 

coul.d  probably  be  eliminated  if  a  weighting  scheme  for  states  with 

28 

eigenvalues  near  the  Fermi  energy,   similar  to  that  used  by  Hattox, 

were  used  for  locating  the  Fermi  energ^y  in  the  present  calculations. 


C.  2    Total  Energy  and  Cohesive  Energy  as  a 
Function  of  a  for  Copper  (A=Ag) 

The  total  energy  and  cohesive  energy  for  the  six  calculations, 
listed  as  pi-esent  calculations  in  Table  C.  1,   are  given  in  Table  C.  2. 
Also  included  in  the  table  are  the  relative  differences  betv/een  the 

calctilated  values  of  cohesive  energy  and  the  experimentally  determined 

34 
value  of  0.  257  Ry.         The  total  energy  as  a  function  of  a  is  plotted  in 

Figs.   C.  5  and  C.  6  for  these  six  solid  calculations  and  for  six  atomic 

calculations  respectively.     Figure  C.  7  is  a  plot  of  the  cohesive  energy 

as  a  function  of  d  for  the  same  solid  calculations. 

As  was  pointed  out  in  Sec.    3.  1.   the  total  energy  for  both  solid 

and  atomic  calculations  is  very  nearly  a  linear  function  of  a  .     However, 


-U' 


TABLE  C.  2 


Total  energy  and  cohesive  energy  as  a  function  of  a  for 
copper  (A=A  ). 


a 

Meta.1 
(Rv) 

Atom 
(Rv) 

Cohesive  Eneroy 
(Rv) 

Relative  difference] 
from  experimental 
cohesive  energy.  ='■  ! 

0.66667 

-3270.7237 

-3270.4534 

0.2703 

0.05 

0. 70000 

-3276.9131 

-3276.6333 

0.2798 

0.09 

0.70635 

-3278.0937 

-3277. 8122 

0.2814 

0.09 

0.72000 

-3280. 6330 

-3280.3480 

0.2850 

0.  11 

0.72250 

-3281.0982 

-3280.8126 

0.2857 

0.11 

0.77000 

-3289.9519 

-3289.6552 

0.2967 

0.  15 

*  Experimental  value  of  cohesive  energy  is  0.  257  Ry. 


34 


Figure  C.  5.     Total  energy  as  a  fiuiction  of  the  parameter  a.,   for 
metallic  copper    (A.=A  ). 
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Figure  C,6.     Total  energy  as  a  function  of  the  parameter  a,   for 
atomic  copper. 
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Figure  C.  7.     The  cohesive  energy  of  copper  as  a  function  of  the 
parameter  a.  of  the  Xq;  3Tiethod  (A.=A  ). 
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A  =  6. 809a. u. 
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dependence  of  total  energy  on  a  cannot  be  carried  too  far.     Even  so,   the 
cohesive  energy  is  a  smooth  function,    increasing  with  increasing  a. 
If  the  cohesive  energy  curve  given  in  Fig.   C.  7  is  extended  to 
intersect  the  experimental  value  of  0.  257  Ry,   the  intersection  will  occur 
at  a  x^alue  of  a  somewhat  below  0.  64.     This  value  of  a  would  give  Fermi 
surface  "necks"  too  large  in  radius  and  a  pressure  at  the  experimentally 
determined  lattice  spacing  of  well  over  100  Kb,  both  of  which  wo\ii.d  be 
in  very  poor  agreement  with  experimental  results.     This  indicates  that 
the  Xo;  method  (using  a  single  a,  as  in  the  present  calculations)  does 
not  simultaneously  give  cohesive  energy,   pressure,   and  a  Fermi 
surface  that  are  in  the  best  possible  agreement  with  experiment. 
However,   the  calculated  values  of  cohesive  energy  differ  from  the 
experimental  value  by  less  than  16%  throughout  this  range  of  a  .     This 
is  a  small  difference  when  we  consider  that  the  cohesive  energies  given 
in  Table  C.  2  are  the  difference  between  two  calculated  energies,   each 
of  the  order  of  3200  Ry.     It  is  quite  possible,   that  increased  precision  in 
the  present  calculations  could  eliminate  the  large  part  of  the  differences 
betv/een  the  calculated  and  experimentally  determined  cohesive  energy. 
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